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CHAPTER I 


Introduction 


The general initial value problem in the theory of ordinary 
differential equations may be defined as the problem of determining the 
values of the m dependent variables y^"^\ j=0, . . . ,m-1 corresponding to 
some desired value of the independent variable, t , subject to 1 ) the 
ordinary differential equation 



(m) 

y 


f(t,y,y 


( 1 ) 


(m-1), 

.y ) 


( 1 . 1 ) 


and 2) a specified set of initial conditions t , y (t ), y^^\t 

o o o o 

y ® dependent variables may be p-dimensional vectors, 

and hence, (1.1) may be a p-dimensional vector differential equation. 
The state vector of the initial value problem is defined to consist of 
the m dependent variables; whereas, (1.1) is referred to as a dif- 
ferential equation of order m . The right-hand side of (1.1), 
f (t,y ,y^ ^ \ . . . ^ ) , is referred to as the function. The initial 

value problem is solved when the state corresponding to some desired 
value of the independent variable is determined. A numerical integra- 
tion method which solves (1.1) directly, i.e., without reducing (1.1) 
to a set of first order differential equations, is referred to as a 
Class m integration method. 



The initial value problem frequently ocrurs in science and 
engineering. One particular application of the initial value problem 
which is of interest is the satellite orbit problem. The ordinary dif- 
ferential equations which describe the motion of a satellite acted upon 
by gravitational and nongravitational forces are given by Newton’s 
Second Law 


d^ -(2) ur (1) 

dt r 


f(r ^ ,t) 


( 1 . 2 ) 


where p = gravitational constant (defined as the product of 
the universal constant of gravitation, G , and the 
mass, M , of the primary body) the primary body 
r = the position vector of the satellite 

r i||7|| 

= the velocity vector of the satellite 
g = vector of gravitational forces 

n = vector of nongravitational forces 


Equation (1.2) is a second-order, nonlinear, ordinary differential 
equation where the function g(r,t) is a smooth, periodic function 
representing the forces of gravity acting upon the satellite. The 
nongravitational contribution may be discontinuous, e.g., entry and 
exit of the satellite into a planet's shadow will suddenly affect the 
solar radiation pressure force. 
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There are many numerical techniques available for solving an 
initial value, ordinary differential equation, and each numerical 
^ method has a limiting degree of accuracy associated with it. The 

! selection of a particular numerical method is subject to the accuracy 

and cost of using the numerical method, where the cost is usually meas- 
ured in terms of computer time. 


t 

s 

I 

If 

’f 


K 
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The numerical integration techniques may be divided into two 
categories: single-step methods and multistep methods. Single-step 
methods require only the value of the state at one value of the 
independent variable in order to advance the solution to another value 
of the independent variable. Multistep methods generally require 
values of the state or values of the function at more than one value of 
the independent variable in order to advance the solution. The values 
of the dependent variable or function, f , that are used to advance the 
solution are determined at values of the independent variable which are 
referred to as the nodes. The stepsizes used to determine the spacing 
between the nodes comprise the mesh sequence. Both techniques require 
the evaluation of the function, f , with values of the dependent vari- 
ables which correspond to various values of the independent variable. 

Equation (1.2) may be reduced from a set of second-order 
ordinary differential equations to a set of first-order equations, thus 
allowing the option of integrating (1.2) with a Class I or a Class II 
method. However. Krogh (1970) and Solis (1975) indicate that (1.2) may 
be integrated more ef.*‘iciently with a Class II method than with a Class 
I method. Since f is a well-behaved, periodic function when 



gravitational forces are dominant, multistep methods are generally more 
efficient than single-step methods from an accuracy versus number of 
function evaluations point of view. Lambert (1973) discusses some of 
the advantages and disadvantages of using a multistep method instead of 
a single-step method. 

Unlike fixed-order/fixed-mesh multistep integrators, 
variable-order/variable-mesh methods estimate the local truncation 
error at each node of the integration in order to satisfy specified 
tolerances. Thus, if a variable-step/variable-mesh integrator uses 
approximately a constant stepsize and order, then a fixed-order/fixed- 
mesh integrator should require fewer computations to integrate the same 
problem with comparable accuracy. For the satellite problem as 
described by (1.2), the fixed-order/fixed-mesh integrators may have a 
significant advantage over variable-mesh/variable-order integrators in 
the amount of computations required, particularly if the orbital eccen- 
tricity is small. 

The most common algorithm for solving the satellite problem 
with a multistep integrator is the PECE algorithm. In the PECE algo- 
rithm, the solution at tj is extrapolated, or predicted (P), forward 
to tj^^ . The predicted solution at tj^^ is used to evaluate (E) 
the function, f . Using this evaluation, the extrapolated solution is 
corrected (C), and a second evaluation (E) of f is made with the 
corrected solution at . Other multistep algorithms include the 
PE(CE)" and P(EO" algorithms where n indicates the number of times 
the steps in parenthesis are applied. However, Krogh (1970) notes that 



for satellite orbits with small eccentricities. 


the Class II/PECE 


methods are more efficient than Class II/PEC methods or Class I/PECE 
methods. 

It is common in the literature to refer to all Class II, mul- 
tistep formulations as Cowell methods and to refer to Class I methods 
as Adams methods. However, at least three distinct Class II methods 
and two distinct Class I methods are available. To avoid any ambigui- 
ties, the terminology of Hersman (1965) will be adopted in which the 
Class II methods are referred to as the general, second-sum, and 
Stormer-Cowell formulations, while the Class I methods are referred to 
as the Adams-Bashforth-Moulton and the first-sum formulations. 

The purpose of this report is to examine the use of two Class 
Il/fixed-mesh/fixed-order /multistep integration packages of the PECE 
type for the numerical integration of (1.2). These two methods are 
referred to as the general and the second-sum formulations. Chapter II 
discusses the derivation of the basic equations which characterize each 
formulation and discusses the role of the basic equations in the PECE 
algorithm. Chapter III discusses possible starting procedures which 
may be used to supply the initial set of values required by the fixed- 
mesh/multistep integrators. In Chapter IV, the results of the general 
and second-siuD integrators are compared to the results of various 
fixed-step and variable-step integrators. 



CHAPTER II 


PECE Algorithms for the General and 
Second-Sum Formulations 

In this chapter, two fixed-mesh/multistep PECE algo, ithmr. ar •• 
developer' for the numerical integration of second-order, initial value, 
ordinary differential equations represented by (1.1). These two formu- 
lations are the general and the second-sum formulations. The basic 
assumptions in developing these fixed-mesh/multistep integrators are: 

1) the value of the state vector (y(t ), y^^^t )) and the value of 

n n 

t are known; 2) the function values f. * f(t,,y,,yi^ , 
n J J J j 

Jsn,n-1 , . . . ,n-i+1 , are known at i distinct nodes where f represents 
the right-hand side of (1.1), and 3) the nodes, tj , associated with 
the values fj are known and satisfy the condition tj * tj_.| ♦ h , 

where Jsn,n-1 n-i>1 and h is a constant ste-'size. The first 

assumption is satisfied by definition of the init. j 1 value problem. 
The second and third are assumed to be satisfied by appropriate values 
supplied from some starting procedure (see Chapter III). 

The development of the PECE algorithm begins by deriving the 
basic equations for each of the formulations. The basic equations 
developed in Sections II. 1 and II. 2 are used to extrapolate or interpo- 
late the state at t from the state at t^ by the proper choice of 

coefficients, if the value of t is between t and t . , , i.e., 

n n— x^i 


6 



if f 


t—t V “t- 

n n-i+i , n 


> - - , then the bajic equations interpolate the state 

, then 


h ' h 

at t , and if the value of t is not between t and t 


n n-i+l 

th' basic equations extrapolate the state at t . Using the extrapola- 
tion and interpolation capabilities of the basic equations, the PECE 
algorithms are determined. A discussion of some modificiations made 
available by the use of back differences comiludes the development of 

the PECE algorithms in this chapter. 


II. 1 General Formulation 

The solution of the second-order ordinary ditferential equa- 
tion given by 




is 


y * y„ ♦ hy^'^ * f f(x,y,y^^^)dx dx^ 

y^^^ * f(x,y,y^^^)dx 


( 2 . 1 ) 


( 2 . 2 ) 


In most applications, the function f cannot be readily integrated by 
analytical means. Thus, f can be replaced by an approximating func- 
tion which represents it to some specified degree of accuracy. The 
derivation of the general formulation algorithm uses the available 
nodes tj and function f(t,y,y^^^) to form a polyncmial P(t) which 
is assumed to yield f. when evaluated at t. . Hence, 
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f(t,y,y<’>) P(t) = 

The polynomial, P(t) , may be written in divided difference form as 


P(t) s f (t— t )g[f,f -] ••• + (t— t )(t— t -) ••• 

n n n n-i n n-l 


(t-t . -)g[f ,f ,,...,f . ,] 

n-i+2 ® n n-1 n-i+1 


^?.3) 


where g[ ] ia the divided difference operator which is defined by 


g[x^,XA,***fX 3 - 

I c n 


«'»i V 

*r*n 


and the degree of the polynomial is (i-1). Imposing the fixed-mesh 
criterion, the divided differences may be written as back differences, 
and the polynomial, P(t) , becomes 


P(t) s f _ ♦ 


(t- 1 ) 


Vf_ ♦ 


n 1! h ”n 


(i-D! " 


« 2 y. (t) 

Jsl J " 


(Z.IJ) 


where 


y^ * 3; yj * 


(t-t . ,) 

(J-DI 


J*2 tStveefi 


and where the back difference operrtort , la defined by 


V® r S f.. and 


k *k 


k-1 


« 2 (-1) f 

a«C 




.J 



-V 
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where | ) are binomial c efficients. 


Using (2.«), it follows that (2.1) and (2.2) become 


y = y, 


♦ + A f' I y.(x)V‘^"''f 


tn o j- 


n <^*1 


(2.5) 


n 


( 2 . 6 ) 


Since only the coefficients, Xj , are a function of the independent 
variable, (2.5) and (2.6) may be rewritten as 


y(t ) = y = y + h^y^^^ + h? 2 a (2.7) 

' n+r 'n+r 'n In I j,r n 


j-1 


' n*r 'n+r 'n I j,r n 


( 2 . 8 ) 


where 


1 s. (shT.)(shT+h)...(shT+(j-2)h) 

^ ! t 4 d, ds, 


o o 


(j-D! h 


j-1 


^ (s,)(s,+h) . . .(sh,+( j-2)h) 

r = ^ 3^1 ““ 

o (j-D! h"^ 


for j s 2,3i...i 


“l,r = i ®1.r * ' "hT ’ 


hr = t -t and t = t ♦ rh. 
I n+r n n+r n 


As defined by Shampine and (k>rdon, basic equations which use i nodes 
are of order i ; thus. Equations (2.7) and (2.8) are defined to be 


order i 
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Equations (2.7) and (2,8) arp thr* equations for solving 
a second-order ordinary differential equation using the general formu- 
lation. Equation (2.8) can be used for solving a first-order ordinary 
differential equation. Krogh (1970) and Sharoplne and Gordon (1975) 
discuss the idea of extending the set of basic equations to solve dif- 
ferential equations of order greater than two. 


Equations (2.7) and (2.8) show that the problem of evaluating 
the integrals of (2.1) and (2.2) has been reduced to a problem of 
evaluating the coefficients, a and B . Sharapine and Gordon (1975) 
present the derivation of the algorithm that is used to evaluate the 
coefficients for a variable-mesh/multistep algorithm. The algorithm 
for calculating fixed-mesh coefficients is a simplification of the 
variable-mesh coefficient algorithm and is presented in Appendix A. 


It should be noted that even though (2.7) and (2.8) are writ- 
ten in terms of the back differences, they may also be written expli- 
citly in terms of the function values as 


n+r 


= y. 




♦ h: 


1 

2 


J.r n+1-j 




•*• h. 


1 

2 6 

Js1 


« 

J.r 


^n+1-J 


(2.9) 

( 2 . 10 ) 


J.r 


(- 1 ) 


J-1 


i 

1 

q=J 



a 


q.r 


where 
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Thus, the basic equations of the general formulation may be written in 
terms of back differences as in (2.7) and (2.8) or in terms of function 
values as in (2.9) and (2.10). 


II. 2 Second-Sum Formulation 


In order to obtain the basic equations for the second-sum 
formulation, certain operators and relationships are required. The 
necessary operators are: 

1 Ic 

1. the back difference operator, V • where V f=Vf-Vf, 

n n n-1 

)c 

2. the shift operator, E , where E f ~ f(t ■•■kh), and 

n n 


3. the differential operator, D , where 


rX d*" ^ (1) dx d^x 

dt*^ dt^ 


( 2 . 11 ) 


As shown by Hildebrand (197^), the above operators satisfy the follow- 
ing relationships: 


e"^° = 1 - V = E"^ 


( 2 . 12 ) 


and 


.-1 


-h 

Ind-^ 


( 2 . 13 ) 




0 


The use of the V operator implies that the stepsize h is a con< 
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stant. 


The derivation of the basic equations begins by making use of 
the three operators above and the relationships from (2.11), (2.12) and 
(2.13). From (2.12), we may write 


'n+s 






'n 'n 


(2.14) 

(2. 15) 


Using equations (2.11) and (2.13), we find that 


= D-"y<2) 

'n 




-h 


2 

In(l-V) 


(1) . p-1 (2) . r -h 

^n I ln( 1-^ In 

By combining (2.14) and (2.15) with (2.16) and (2.17), 
following relationships 


(2.16) 

(2.17) 


we obtain the 



( 2 . 18 ) 


(2.19) 


By expanding the terms in the brackets of (2.18) and (2.19) in an 
Infinite series, i.e.. 


(1-7)~* 

(ln(1-7))^ 


2 

J=-i 



( 2 . 20 ) 



oo 


(1-V) 


-5 


ind'::^ °j,s 


V«Jf 


(2.2V 


it follows that (2.18) and (2.19) may be rewritten as 


n+s 


= h 


j=-2 




y 


(1) 

n+s 


h I 
J=-1 


b 


J.s 



( 2 . 22 ) 

(2.23) 


The coefficients a. and b. are defined in terms of s by 

J.s J.s 

expanding the left-hand sides of (2.20) and (2.21) in an infinite 
series in V and comparing the coefficients of the different powers of 
the V operator. Derivation of the recursive algorithm used to calcu- 
late the coefficients is cumbersome and is discussed by Spier (1971) 
and Velez and Maury (1970). The algorithm used to generate the coeffi- 
cients for the fixed-step second-sum formulation is summarized in 
Appendix B. It should be noted that the coefficients a ~ ® i 

and b . are such that 
-1,s 


a 


-2.S 



1 


a 


-1,3 


S-1 


”jr all s . Thus, (2.22) and (2.23) may be rewritten as 
2 


y s h 
'n+s 




V-2f - (1-s) a, , vJf 

n 11 J~0 ^ J 

I V’V ♦ V b, . vJf 
L " jso . 


(2.24) 

(2.25) 





14 


where the sucimation terms have been truncated to include the first i 

terms. Equations (2:.2'0 and (2.35) are said to be of order i . In 

-1 -2 

(2.24) and (2,25). the terms V f^ and V f^ are referred to as the 
first and second terms, respectively. The first and second terms 

satisfy the following relationships. 


v'V 

= vV 

+ f 

(2.26) 

n 

n-1 

n 


V-2f 



(2.27) 

n 

n 



from the definition of the back difference operator. It should be 
noted that only (2.25) is required to solve a first-order ordinary dif- 
ferential equation. 


By using function values instead of back difference, (2.24) 
and (2.25) may be rewritten as 


y 


n+s 


y 


( 1 ) 

n+s 


1 V"\ - (l-s)V’f a* f , I 

♦ % V, Vj ] 



(2.28) 

(2.29) 


where 
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II. 3 Relati onships Between the General and Second-Sum Coefficients 


To illustrate the relationship between the coefficients of 
the general formulation and second-sum formulation coefficients, each 
set of basic equations for extrapolating the solution one step is 
presented and compared. To predict the solution at t^^^^ from t^ , 
the general formulation basic equations are 


'n+1 


y„ ♦ h,<’> * 2 c vJ-'f 

n n ^ j t ' ^ 


'n+1 


y- ' + h 2 e, , 


'n 


j=1 


j.1 


(a. 30) 
(2.31) 


The form of the basic equations for the second-sum formulation that 
predict the solution at t^^^ from t^ are 


y , s h^ f V"^f ♦2a.,, V’^"V 1 

n+1 I n jti "J 

y^’J = h [ V“V ♦ 2 b, , , 1 

n+1 I n jsi nj 

and that interpolate the solution at t from t are 

n n 

y = h^ [ V"^f - V"V ♦ 2 a, , 1 

'n I n n J-1,0 nJ 

y^^^ = h f V"V ♦ 2 b, , n V'^”^f 1 . 

" n J-1,0 n J 


( 2 . 32 ) 


(2.33) 


(2.3*0 


(2.35) 


From (2.31), it follows that 
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= h 5 e . 

n+1 j,i n 


Using (2.33) and (2.35), we find 




and comparing this equation to (2.36) t we also find 


V) ' Vi.i " Vi.o ■ • 


(2.36) 


(2.37) 


Hersman (1965) has shown that, in general, . 

Thus, (2.37) becomes 



= b 


J.l 


or 


’j,k ^ ^J,k 


Similarly, by forming the back difference Vy^^^ for each formulation 
and comparing coefficients, we also find that 


“j.1 * ■ Vi.o 


' ' Vl.O 


(2.38) 


Hersman (1965) has also shown that * *k,1 * (2.38) may 

be written as 


' *J.1 ' Vi.o 




and 



oil, = a,, - b. 

J.k J.k j-1.0 


II. *• Development of the PECE Algorithm Equations 


The fixed-mesh, PECE algorithm assumes that the values fj 
and tj , for j=n,n-1,...,n-i-*-1, are knovm, that t^ satisfies the 
condition t. . = t ♦ h , and that if t, y and are known, then 
f = f(t,y,y^^^) can be calculated. The PECE algorithm is based upon 
the extraf>olation property of the basic equations. The basic equations 
for the general formulation are 


, . y . h,<’> . i a, vJ-’f 

P n 'n j,s n 


= y ♦ hy^’^ + h^ 2 f 4 i 

n ^n j,s n-j+1 


( 2 . 39 ) 


y<’> = yi’>.h i 6, 




( 2 . 40 ) 


and for the second-sun formulation are 


yp-h 




( 2 . 41 ) 
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= h 






( 2 .*» 2 ) 


where s = p-n. 

The prediction step of the PECE algorithm is an extrapolation 

of the solution from t to t , . By setting s=1 in (2.39) through 

n n-f 1 

(2.42), the predicted solution at t^^^^ , i.e., ^n+1 ’ 

given by the general form of the basic equations as 


'n+1 


. (1) V.2 i 

= Yn ^ hyn ♦ h ^2^ ^ f„ 


.( 1 ) .2 


n^1 


y'” . h 2 6 7l-’f 

" J = 1 " 


( 1 ) 


• 'n * " J/j.t Vj*' 


and by the second-sum form of the basic equations as 


p , = h 


♦ 2 a, , , V« -f 


vJ-' 


j, 

* jf , V '.' ] 


- h 

Pn^l * ^ 


1 


(2.43) 


(2.44) 


(2.45) 


(2.46) 
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Usinp the values for t . , , 

n+i n*#*! 


and 


Vi • 


the first 


( 1 ) 


evaluation of the PECE algorithm is given by fj^^^ = ’^n+1 


) . 


A new set of back differences, 

Vi-j Vi • 

defined by 


V'^f* , , or function evaluations, 
n*M 

The new set of back differences are 


V°f* , = fP 
n+1 n+1 


V^f* , = , 


n+1 


‘n+1 


J=1 1 


and the new set of function evaluations are defined 


f = fP 
‘n+1 n+1 


fj s f^, j=n,n-1,...,n-i+1 . 

With the inclusion of fP^^ , there are i+1 back differences or func- 
tion evaluations available to approximate the function f . The new 
polynomial approximation of f now spans from ^n+1 

instead of t. . - to t, . 

n-14-1 n 


The predicted solution, p^^^ and p^|j , is corrected by 
extrapolating the solution at t^ , i.e., y^ and y^^^ , to t^^^ by 
use of the new polynomial order i . The basic equations give the 
corrected solution as 


Vi 


. (1) k 2 oJ-V* 

^n ^n J.o n+1 


Ml 2 • • 

‘ '”'1. * " “j.« Vz-J 


s 


(2.M7) 
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( 1 ) 


(1) , 

♦ h 2 


'■'♦h 1 r 1 

n j.o n+1 


( 1 ) • • 

' ^ ^J.o ^n^2-j 


for the general formulation and 


2 1-2 -1 * • 
y„ 1 = h V f - V f ♦ 2 a, , Vj-1f , 

n+1 I n n j-ltO n+1 


2 -1 a • 

^n " ^n ”*■ ®J-1,o *^n+2-j 


I 

...[ 

- h [ v"V ♦ 2 b v^"V* 1 

'n +1 - I '' n jfi J- 1.0 n +1 J 

r -1 • • 1 

[' ^ I Vl.o ^n.2-J J 


= h 


for the second-sum formulation. 

Tne second evaluation of the PECE algorithm is made by 
^.he corrected solution, i.e., 

*^n+1 * ^n*r * 

The back differences are updated by the fornulas 


s ix1 i 

'^*^n+l Vi *^n * ^ ^ ‘ 


and the first and second sums are updated by 


, * '/“V ♦ f , 

n^1 n n+1 


(2.4b) 


(2.49) 


(2.50) 


using 
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v-^f 


n+1 


= 




n+1 


II. 5 A Modification to the PECE Algorithm 

Tlie use of back differences instead of function evaluations 
in the PECE algorithm allows for a simplification of the correction 
formulas. The simplification will be demonstrated here as it applies 
to the general formulation, but it may also be applied to the second- 
sum formulation. The process of propagating the back differences is 
also discussed in this section in order to complement the use of the 
back difference form of the PECE algorithm. 


Equation (2.3) states the approximation of the function f 
as a polynomial in t and using i nodes and i function evalua- 
tions. In divided difference notation, the polynomial approximation of 
f was given as 


P(t) = f ♦ (t-t )gtf ,f ,] ♦ ••• 
n n n n-i 


. <t-t„) - CJ.51) 


where g[ ] is the divided difference operator. The general formula- 
tion basic equation was derived by writing (2.51) with back differences 

as 


I 

<2.52) 

jsi j n 

and integrating the coefficients, YAt) , from t to t . 

j n 
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Equation (2.51) rf rr^'r ? the polynomial used in the predic- 
tion step of the PECE algorithr.. The polynomial used in the correctior 
step of the algorithm uses the i nodes and i functicn evaluations 
used in the predicting step and an additional function evaulation and 

node, f^ , and t , . The polynomial that incorporates the i-*-1 
n^l n+’i 

nodes and i+1 function evaluations may be written as 

P*(t) = r ♦ 1^ 

n n n n - 1 


♦ (t-t^)...(t-t^_^^2^*^^n'^n-l** 

= P(t) ♦ (t-t )...(t-t^ 

n n-i>i n 


• • t ^ . « 1 

n-i>1 


f f^ 1 

’ n-i+r n>r 


f 1 

n-i+1’ n+r 


(2.53) 


Equation (2.53) represents the polynomial used in the correction step 
of the PECE algorithm. The corrected solution is given by the calculus 
solution 


'n+1 


^ (1) 

^n " “^n " / / 

^n ° 


^n+l 


, P*(x)dx 

t 

n 


P*(x)dx dx^ 


(2.5«) 

(2.55) 


Using (2.53) with (2.5*0 and (2.55), we find 
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y 

«. ■ 


( 1 ) ,^ n+1 ,*1 

- J j P(x)dx dx^ 

t o 

n 


'n+1 - ^ "'n 


» /"*’ /’ dx, 

t o 
n 


= p„»i * >'1.1.1 ’'Cl 


(2.56) 


. 0 ) 


'n+1 'n 


y; ' + P(x)dx + 


n+1 


n 


n 


X pO; . h e. , , v^fP , 

n+1 i+1,1 n+1 


(2.57) 


where 


vV , 

n+1 


Cl 


- 2 V‘J~V, 

j=i 


(2.58) 


Equations (2.56) and (2.57) represent the correction formulas that 
would be used in place of (2.47) and (2.48) for the general formula- 
tion. Similarly, for the second-sum formulation, (2.49) and (2.50) 
could be replaced by 


y 

y 


n+1 

( 1 ) 

n+1 


**n+1 ^ ®i+1,1 ^ Cl 


Cl * ^ **1+1,1 ^Cl 


(2.59) 

(2.60) 


One advantage of (2.56) and (2.57) is that, instead of 21, 
only i+1 coefficients are needed to predict and correct the solutioh. 
Thus, (1-1) fewer multiplications per step are required. The addi- 
tional computations required to form and for propagating the 
back differences from t^ to t^^^ are a possible disadvantage to 
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(2.56) and (2. 57) However, the propagation of the back differences nay 
take advantage of the internediote calculations resulting from the cal- 
culation of . By retaining the intermediate sums, 

^ 1-1 

S. = I f k=1 i 

j=k " 

resulting from the evaluation of (2.58), the back differences at 
are then calculated from 


V^f , = S r f , - S, 
n^l o n+1 1 


''“■''■n.l = \ \ ‘ • 


The modified PECE algorithm is given in Appendix C for the general for- 
mulation and in Appendix D for the second-sum formulation. 



II. 6 Notes and Comments 

The algorithm discussed in this chapter is for a PECE mul- 
tistep integrator with a predictor of order i and a corrector of 
order i>1 . The PECE algorithm is one of a family of algorithms in 
which corrector formulas and the function evaluation using the 
corrected solution may be applied any number of times. This family of 
algorithms is represented by PE(CE)" and P(EC)". The choice of algo- 
rithm to be used to solve a given problem is discussed by Shampine and 
Gordon (1975) and Krogh (1970). 


The order of the correcting formula is generally chosen to be 
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of the same order or one higher order than the predicting for'"ula. The 
decision of which order to use for the corrector is dependent or. the 
differential equation to be solved and on the order of the predictor. 
Considerations for selection of the corrector order are discussed by 
Shampin® and Gordon (1975) and Krogh (1970), 

The fixed-mesh/raultistep algorithm advances the solution in 
intervals of the stepsize, h , of the independent variable t . How- 
ever, the solution may be desired for some value of t , e.g., tj , 

such that t. is between two nodes. The solution at t. is found by 
J J 

advancing the solution until tj is bounded by two of the nodes. If 

the nodes range from t to t . , such that 

n n— 1+1 

It - t . , 1 > It, - t 1 , the solution at t is interpolated by 
n n-i+i — j n 

using r = s = j-n in the basic equations. 

The problem of obtaining the Initial set of function evalua- 
tions for the PECE algorithm is discussed in Chapter III. The general 
and second-sum algorithms discussed in this chapter are compared to the 
other numerical integration packages in Chapter IV. 



CHAPTER III 


Starting Methods 

As discussed in Section II. 4, a multistep numerical integra- 
tion algorithm of order i assumes that i function evaluations and 
i nodes are known. The initial values of the nodes and function 
evaluations are found by application of an appropriate starting pro- 
cedure. A starting procedure has only the initial conditions, , 

y and y^^^ . and the differential equation, y^^^ = f(t,y,y^^^) , 
o o 

available to calculate the values of the nodes and to evaluate the 
function. Some of the various starting algorithms, which can be used 
with roultistep integrators, include 1) the bootstrap method, 2) the 
iterative method and 3) the use of a single-step integrator. This 
chapter begins by discussing the starting procedures for the Class II 
fixed-mesh/ fixed-order multi step integration packages using the general 
and second-sum algorithms described in the previous chapter. The 
chapter concludes by detailing a proposed starting procedure for the 
Class Il/fixed-mesh/fixed-order/multistep integration packages which 
are the emphasis of this report. 
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III. 1 Summary of S tart Ing A lgorithms 

The bootstrap starter begins with a first-order ba..ic equa- 
tion and increases the order as the solution is advanced one node at a 
time. By using the initial conditions, the first function evaluaL' >r. 
is made, f = f(t ,y . The solution is extended to t. = t + h 

O O O O i 

by using a first-order basic equation, e.g., i=1 in (2.39) and (2.40) 
or in (2.41) and (2.42). With the predicted solution at t^ , another 
function evaluation is made, = f(t^,p,,p^^ ^) . The correction form 
of the basic equation is employed and an additional function evaluation 
is made with the corrected solution at t^ to complete the PECE algo- 
rithm at t^ . Now, with two function evaluations available, the solu- 
tion is advanced to t 2 by application of the PECE algorithm for 
i=2 . In this manner of bootstrapping, the solutior is advanced until 
the required number of function evaluations are known. This starting 
procedure is generally performed with variable-step/variable-order mul- 
tistep integrators in which the local truncation error is estimated 
after each step in order to determine if the step is acceptable and to 
determine the next stepsize or order to use. 

A single-step numerical integrator, e.g., a Runge-Kutta 
method, may be used to generate the required function evaluations. By 
using the single-step integrator to advance the solution from one node 
to the next, the solution at the nodes is obtained and the function 
evaluations are calculated and stored. 
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An iterative starter assumes that a first approximation of 
the required function evaluations has been made by means of a bootstrap 
method, by use of a single-step integrator, by a Taylor series expan- 
sion or by any other suitable procedure. Once the first approximations 
have been provided, the basic equation is used to interpolate the solu- 
tion at each of the nodes, and a second set of function evaluations is 
calculated. The second set of function evaluations is used with the 
basic equation to interpolate the solutions at the nodes again and a 
third set of function evaluations is calculated. The iterative method 
proceeds in this manner until some termination criteria is satisfied. 
The criteria may be a certain number of iterations, or it may depend 
upon the difference between the solutions of two consecutive itera- 
tions. Once a criterion is satisfied, the required function evalua- 
tions are known. 

III. 2 The Central Iterative Start Ing Algorithm 

Two important criteria were considered in selecting a start- 
ing procedure for the general and second-sum packages of this report. 
The first criteria fo.' '.he starting procedure was to generate a set of 
function evaluations and nodes wit'' '«''ect to the initial conditions 
which can be used to interpolate the initial conditions exactly. 
Secondly, a starting procedure should be consistent with the interpola- 
tion scheme of the integrator. The interpolation scheme chosen for use 
in the fixed-step/fixed-order integrators of this report requires that 
the solution be advanced far enough that the point at which the 
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interpolation is to be performed is approximately midway between the 
extreme nodes. This scheme helps to reduce discontinuities in the 
interpolation when the solution is advanced from one node to the next 
node. The starting procedure used in the general and second-sum 
integration packages described in this report is a central, iterative 
starting procedure. 


Backward or forward iterative starting procedures are such 
that the initial set of nodes span the interval t to t, , or the 
interval t^ to i • However, *> central starting procedure selects 
the end nodes such that the initial conditions lie at a node approxi- 
mately midway between the end nodes. With the initial conditions at 
t^ , the end nodes will be at tj and t^ , where q = integer (^) 
and J = q - i + 1 . Thus, the central starting method will have 

advanced the solution to t and will return the function evaluations 

Q 

at the nodes t^^; k=q,q-1 , . . . ,q-i+l . 


As mentioned above, one important requirement for the start- 
ing procedure is that the solution at t^ and the function evaluations 
must be consistent with the initial conditions. For the general formu- 
lation using back differences, this imposes the condition that 


= y^ - qhy^^^ ♦ q^h^ 5 


i 

5 

J=1 


o 'q 

( 1 ) ( 1 ) 

yo = y, - qh 


. v-'V 

J.-q q 


5 B, _ 


( 3 . 1 ) 

( 3 . 2 ) 


The solution at any of the nodes used in the starting procedure is 
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given by 

(3.3) 

(3.4) 

for k=q,q-l , . . . ,q-i+1 . 




y = y + (k-q)h + (k-q)^h^ 2 ol . 

k 'q q J.k-q q 


( 1 ) ( 1 ) 


= y^ + (k-q)h 2 B 

j -1 


j.k-q 


vJ-^ 


The interpolation formulas for the central iterative starter may be 
found by using (3.1). (3.2), (3.3) and (3.4) to obtain 


Vi. = 


Yo khy^^^ 


+ 2 ( 
j=1 \ 


(q-k)^a 


. . - Q ot, kq 6. 

J.k-q ^ j,-q ^ j 




(3.5) 




(3.6) 


k=q,q-1 .... ,q-i+1 . 

Equations (3.5) and (3.6) use the initial conditions at t^ and the 
back differences at t^ to interpolate the solution at the nodes. 
Equations (3.5) and (3.6) will always satisfy the initial conditions. 


The initial values of the function evaluations are generated 
by using the Taylor series expansions 

t, r t + kh 
J o 
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y, = y„ khy^'^ + k^h^f 
J o o c 

( 1 ) ( 1 ) , 

»k -■ >'o * <^0 






to approximate the solution at the nodes in order to evaluate the func- 
tion f(t,y,y^^^) . Using the initial values of the function evalua- 
tions, the iterative procedure uses (3.5) and (3.6) to find the final, 
required set of function evaluations. 


Since the function f is assumed to be smooth and continu- 
ous, the solutions at the nodes during the iterative procedure are 
assumed to converge upon a trajectory. The final trajectory, the accu- 
racy of the trajectory and the number of Iterations required to con- 
verge on the trajectory are a function of the stepsize h , the order 
i and the function f . The starting algorithm employed by the gen- 
eral and second-sum packages measures the convergence of the solution 
by calculating the relative norm of the difference between two succes- 
sive solutions at t . If z, is defined to be a vector composed of 

Q J 

the dependent variables y and y^'^ at t * t^ on the itera- 
tion, then the relative norm, u , of the difference between two suc- 
cessive iterations at t is defined to be 


2n, 


j? 




where, for k s 1 , . . . , 
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z.(k) - 2 , (k) 

z(k) = nrf-^ ^ 0 

ZjCK; J 

= Zj_^ (k) if = 0 , 

and n is defined to be the number of elements in each of the vectors 

y and • The solution at t^ is assumed to contain the greatest 

error and to be the slowest converging solution at any of the nodes 

because t is generally the furtherest node from t . 

q o 

Equations (3>5) and (3.6) are the equations used in a start- 
ing procedure with the initial conditions at t . This notation does 

o 

not lend itself to simple application in a computer program. Appen- 
dices E and F give the general and second-sum equations equivalent to 

(3.5) and (3.6), but where the initial conditions are at t and the 

n 

end nodes are at t and t , , where t is between t and 

n n-i+1 ffl n 

t , , . It should be noted that the use of (3.5) and (3.6) depends 
n— I 

upon the calculation of a special set of coefficients; however, the 
calculation of the coefficients require little effort and are only 
required for the starting procedures. 


» 



CHAPTER IV 


Comparisons and Results 

The evaluation of the general and second-sum fixed- 
mesh/fixed-order Class II multistep integrators was carried out in two 
phases. The first phase compared the performances of the general and 
the second-sum integrators to the performances of four Class I integra- 
tors and two Class II Integrators in solving a selected group of dif- 
ferential equations with periodic solutions. The second phase compared 
the transverse errors of these integrators for two typical satellite 
problems with the following gravitational force models. Da spherical 
earth modeled as a point mass and 2) a non-spherical earth modeled by 
an eleventh degree and order spherical harmonic geopotential. The 
results of the first and second phases are discussed in Sections IV. 1 
and IV. 2, respectively. 

Three forms of the flxed-mesh/fixed-order Class II multistep 
integrators were adopted in this investigation. The general formula- 
tion of the Class II integrator package, referred to as KSGFS, used the 
back difference form of ''he PECE algorithm. TWo forms of the second- 
sum formulation were used. The first second-sum package, referred to 
as SSFSBD, used back differences while the other second-sum package, 
referred to as SSFSFE, used the function evaluation form of the PECE 
algorithm. 
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The reruli’? of the peneral and the second-sum Class II 
integrators are compared with three documented integration packages and 
two undocumented Integration packages. The documented software pack- 
ages were ODE by Shampine and Gordon (1975), KROGH by F. T. Krogh 
(1970) and RK(7)8 as derived by Fehlberg (1968) and applied by McKenzie 
and Sepehnoori (1978). It should be noted that KROGH may be used as a 
Class I, Class II or a combination of both Class I and Class II 
integrators. To distinguish between the first two of these modes, the 
Class I ^orm and the Class II form of KROGH will be referred to as 
KROGH 1 and KR0GH2, respectively. The undocumented packages were ABFS 
and RKN7(8). ABFS is the Class I equivalent of KSGFS. RKN7(8) is a 
variable-step Runge-Kutta-Nystrom integrator used for the solution of a 
general ordinary differential equation of Class II. The coefficients 
for RKN7(8) were derived and published by Fehlberg (1974). 

All computer work was carried out on a CDC Cyber 170/750 com- 
puter using the FTN compiler at the University of Texas at Austin under 
the UT2D ooerating system. The Cyber 170 computer utilizes a 60-bit 
word; 12 bits are used for the sign and exponent, and 46 bits are used 
for the mantissa which results in about 14 decimal digits of accuracy. 

IV. 1 Integrator C< -Tipariaons for Differential Equations 
with Periodic Solutions 

The first phase of the evaluation of KSGFS, SSFSBD and SSFSFE 
consisted of comparing their results with those from ODE, KR0GH1, 
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KR0GH2, ABFS, RFC(7)8 and RKN7(8) for the solution of five sets of dif- 
ferential equations that have periodic solutions. The differential 
equations and results are discussed In Sections II. 1.1 through II. 1.5. 


The comparisons of the integrators for this phase was accom- 
plished by using the software package COMPAR-* COMPAR allows the user 
to plot the efficiency curves for each integrator or combinations of 
integrators. The efficiency curves are defined as the endpoint error 
versus the number of function evaluations required, endpeint error 
versus the amount of central processor time required, the maximum glo- 
bal error versus the number of function evaluations and the maximum 
global error versus the amount of central processor time. The global 
error at each step is defined as 


GE = Maximum 


x^(i)-x_(i) 
C i 

b 


J 


i=1 


.n 


whers n is the number of elements in the state vector x , x is the 

c 

state vector calculated by the integrator, Xi^ is the reference solu- 
tion and b is the maximum (x,p(i),1.). The endpoint error is the glo- 
bal error at the final point of the integration interval. By varying 
the tolerances, the efficiency curves of the variable-step integrators 
could be determined. To obtain the efficiency curves for ABFS, KSGFS, 
.SSFSBO and SSPSFE, several computer runs were made to determine the 
optimum order for each integrator on each problem. By using the 

* CONPAR was developed at the Department of Aerospace Engineering and 
Engineering Mechanics, The University of Texas at Austin, by 
Richard McKenzie. COMPAR is used to compare one or more 
integrators, each at one or more tolerances, for solving a set of 
differential equations. 
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optimum orders and varying the stepsize, the efficiency curves for the 
fixed-mesh/fixed-order integrators can be found. It should be noted 
that the starting procedures for ABFS, KSGFS, SSFSBD and SSFSFE we< 
allowed to converge the relative error norm (see Chapter III) as far as 
possible. Generally, the relative error norm was reduced to zero. 

Table iV.1.1 through IV. 1.5 are 'ummaries of the results from 
COMPAR for each of the five problems discussed in Sections IV.1.1 
through IV. 1.5. The first two columns of these tables are labelled 
either ABSERR and RELERR or ITER and BORDER. For the variable-step 
integrators, ABSERR and RELERR are the absolute and relative error 
tolerances to be used by each integrator. For the fixed-mesh multistep 
integrators, ITER represents the maximum number of iterations to be 
used in the starting procedure, and BORDER represents the order of the 
integration formula to be used. The next three columns in the tables 
are labelled BFE, BSTPA and BSTPR and are the tot 1 number of function 
evaluations made, the total number of accepted steps by each integrator 
and 'he total number of steps rejected by RK(7)8 or RKN7(8). The two 
columns labelled CP-TIME and OVHD indicate the amount of central pro- 
cessor time used by each integrator for each set of ABSERR and RELERR 
tolerances. CP-TIHE Is the total amount of central ^<'Ocessor time used 
by the integrator and the derivative evaluation routines, while OVHD is 
tne overhead or the CP-TIHE excluding the time spent in performing 
derivative evaluations. The remaining columns are self-explanatory. 


The problems discussed in Sections IV. 1.2 through IV. 1.5 were 



suggested by Shampine and Gordon (1975). These four problems have non- 
linear differential equations and ’’give a good indication of how the 
codes perform on problems of celestial mechanics."** The harmonic 
oscillator problem of Section IV. 1.1 has the same analytic solution as 
the circular two-body problem of Section IV. 1.2 out has a set of linear 
differential equations. The elliptic two-body problem of Section 
IV. 1.3 has the same set of differential equations as the circular two- 
body problem but has a different analytical solutio.:. The Euler 
rigid-body problem of Section IV. 1.4 has a Jacobian elliptic function 
as the solution. The restricted three-body problem of Section IV. 1.5 
does not have an analytical solution and presents a case of rapidly 
changing derivatives near the close approaches. 


IV. 1.1 Harmonic Oscillator Problem 

The harmonic oscillator problem may be modeled by the set of 
first-order differential equations 


dxi 

It'- *2 


dt ' "*1 


or by the second-order differential equation 


d^x. 


dt 


2 ■ “*1 


•* Shampine and Gordon (1975), page 242 



* 18 
With the initi'l i'.-"nHtions x,(0)=0.0 aiid x-,(0) = 1.0 , the analytical 

I C. 

slut ion is ^'ivcn by 

t 

' x^(t) : sin(t) 

X 2 (t) = cos(t) 
which has a period of 2n . 

Several computer runs were made with the integration interval 
ranging from 0.0 to 2^’, 0.0 to 12 ti and 0.0 to 20 ti. By varying the 

stepsize from 0.5 to 0.3 for these three intervals, the optimum orders 
for the fixed-mesh/fixed-order multistep integrators were determined to 
be 8 for ABFS, 15 for KSGFS and 11 for SSFSBD and SSFSFE. 

The efficiency curves for this problem were determined for 
the integration interval of 0.0 to 20rt . The efficiency curves for the 
variable-step integrators RK(7)8, RKN7(8), ODE, KR0GH1 and KR0GH2 are 
shown In Figures IV. 1.1a through IV.I.Id. Figures IV. I.le through 
IV.I.Ih show the efficiency curves for ABFS, KSGFS, SSFSBD and SSFSFE 

relative to those of ODE, KR0GH1 and KR0GH2. Figures IV. 1.11 through 

IV. 1.11 show the efficiency curves of ABFS, KSGFS, SSFSBD and SSFSFE 

relative to those for RK(7)8 and RKN7(8). A summary of these results 

are given in Table IV. 1.1. 




IV. 1.2 Circular Problem of Two Bodies 

The circular problem of two bodies may be modeled in two 
dimensions by the set of first-order differential equations 








Figures IV. 1.1 a and b 

Efficiency Curves for the Harmonic Oscillator Problem 










Figures IV. 1.1 c and d 

Efficiency Curves for the Harmonic Oscillator Problem 






Figures IV. 1.1 e and f 

Efficiency Curves for the Harmonic Oscillator Problem 
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Figures IV. 1.1 g and h 

Efficiency Curves for the Harmonic Oscillator Problem 
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Efficiency Curves for the Harmonic Oscillator Problem 





Efficiency Curves for the Harmonic Oscillator Problem 
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1 .518E-01 

.19100 

1 .OOE-08 

1 .OOE-12 

827 

422 

0 

.00000 

1 .4891-01 

.19100 

l.OOE-10 

1. OOE-12 

1584 

800 

0 

.00000 

7.854E-02 

.38000 

l.OOE-11 

1. OOE-13 

1600 

80V 

0 

.00000 

7.767E-02 

.38500 

1. OOE-12 

1. 001-13 

1632 

83; 

0 

.00000 

7.561E-02 

.40100 

1 .OOE-13 

l.OOB-14 

1639 

834 

0 

.00000 

7.5341-02 

.39100 


OVHD 


MAXIMUM 
END POINT GLOBAL 
ERROR ERROR 


MINIHL'M MAXIMUM 
STEP SIZE STEP SIZE 


.02661 8.287E-03 B.287L-03 2.000E-02 
.05846 2.432E-05 2.432E-05 2.000E-02 
.09351 2.752E-07 2.752E-07 2.000E-02 
.14358 2.911E-09 2.912E-09 2.000E-02 
.20668 3.186E-10 3.187E-10 2.000E-02 
.27610 3.379E-11 3.361E-n 2 .OOOE-02 
.37192 4.262E-12 4.280E-12 2.000E-02 


2.308E*00 
1 .116E^00 
7.044E-01 
4.178E-01 
3.133E-01 
2.461E-01 
1 .938E>01 


HAXIHim 

END POINT GLOBAL MINIMUM MAXIMUM 
OVHO ERROR ERROR STEP SIZE STEP SIZE 


.01596 

.03130 

.05699 

.09669 

.12110 

.18136 

.21962 


7.330E-04 
7.293E-06 
7.444E-08 
7.554E-IO 
7.639E-11 
1 .OOOE-11 
1.166E-12 


7.330E-04 

7.293E-06 

7.444E-08 

7.554E-10 

7.639E-11 

l.OOOE’lt 

1.838E-12 


2 .OOOE-02 
2 .000E>02 
2.000E-02 
2.000E-02 
2.000E-02 
2.000E-02 
2.000E-02 


1 .661E«00 
8.925E-OI 
4.986E>01 
2.800E-01 
2.098E-01 
1.586E-01 
1.262E>01 


OVHD 


MAXIMUM 

END POINT GLOBAL MINIMUM MAXIMUM 
ERROR ERROR STEP SIZE STEP SIZE 


.06526 2.794E-03 
.09786 1.958E-05 
.14946 2.719E-07 
.20225 3.969E-09 
.26239 5.762E-II 
,29127 4.769E-11 
.40898 2.843E-11 


2.784E-03 

1.964E-05 

2.679E-07 

3.949E-09 

7.906E-07 

7.500E-07 

5.534E-07 


2.500E-03 
2.500E-04 
2.500E-05 
2 . 500E-06 
7.906E-07 
2.500E-07 
7.906E-08 


6.4O0E-0I 
2.560C-01 
2.048E-01 
1 .638C-OI 
1 .036E-01 
1.31IE-01 
1.343E-01 


END POINT 
OVHD ERROR 


MAXIMUM 

GLOBAL MINIMUM MAXIMUM 
ERROR STEP SIZE STEP SIZE 


.08709 2.400E-05 
.10390 5.777E-08 
.17384 5.465E-U 
.19013 6.814E-12 
.18026 6.579E-12 
.33338 1.852E-11 
.36617 2.114E’ll 


2.488E-05 1.000E<02 6.400E-01 
5.716E-08 1.250E-03 3.200E-OI 
5.802E-11 7.813E-05 1 .600E-01 
6.736E’12 9.766C-06 1 .600E-01 
6.499E-12 2.441E>06 1 .600E-01 
1.851E-II 1.22IE-06 1 .600C’01 
6.104E-07 3.052E-07 8.000E-02 


OVHD 


MAXIMUM 

END OINT GLOBAL MINIMUM MAXIMUM 
EMOt ERROR STEP SUE STEP SUE 


.08506 

.17354 

.17318 

.34586 

.35051 

.36582 

.35567 


3.466E-04 3.424R-04 1 .OOOE-02 3.200B-01 
5.165E-08 5.134E*08 1.250E-03 1 .600EH)1 
4.096E<08 4.079E>08 7.813E-05 1 .600E>0I 
5.253E-11 5.250E-I1 9.766E-06 1 .600E-0I 
1.994B<1I l.988E-n 2.441E-06 1.6008-01 
2.163E-11 2.156E-11 1 .221E-06 8.000E-02 
2.150E-11 6.104E-07 3.0521-07 8.000E-02 


Table IV. 1.1 

COMPAR Summary of Statistics for the 
Harmonic Oscillator Problem 



INTM'RATION MrTHaO A»FS 


MAX I HUH 


urn 

ORPFR 

NEE 

NSTPA 

NSTPR 

NSTPR/ 

NSTPT 

AVERAGE 
STEP SIZE 

CP 

- TIME 

OVHO 

END POINT 
ERROR 

1. LORAL 
ERROR 

MINIMUM 
STEP SIZE 

MAMV. 

STEP SUE 

2.50E*OI 

8.00E*00 

2609 

1256 

0 

.00000 

5.003E-02 


20500 

.14077 

9./02E-12 

9 .666E-12 

5.000E-02 

2.500E-01 

2.^0F*OI 

8.00E*00 

1 569 

628 

0 

.00000 

1 .OOIE-01 


.11100 

.08149 

4.951E-09 

4.917E-09 

1 .OOOE-01 

5. OOOE-Ol 

2 .50E»0I 

8.00E*00 

1117 

502 

0 

.00000 

I-252E-01 


.08750 

.06342 

3.679E-08 

3 .648E-08 

1 .250E-0I 

6.250E-0I 

2.i0E401 

8.00E*00 

965 

418 

0 

.00000 

1.503E-01 


.07700 

.05620 

1 .893E-07 

1 .873E-07 

1 .500E-01 

7.500E-01 

2.i0E*0l 

8.C0E*00 

773 

314 

0 

.00000 

2.001E-01 


.05650 

.03984 

2.509E-06 

2.477E-06 

2.000E-01 

1 ,OOOE*00 

2.^0E*0l 

8.00E*00 

655 

251 

0 

.00000 

2.503E-01 


.05600 

.04188 

1 .854E-05 

I .025E-O5 

2 .500E-01 

1 ,250E*00 

2.^0E*0I 

8. 00 r* 00 

595 

209 

0 

.00000 

3.006E-01 


.04800 

.03518 

2.888E-02 

2.999E-02 

3.000E-0I 

1 .500E*00 

INTLCRATIOM METHOD 

KSCFS 






















MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

IT^P 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP 

- TIME 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SUE 

2.50C*0I 

1 .5OE*0l 

2633 

1256 

0 

.00000 

5.003E-02 


.23800 

. 18125 

6.762E-14 

2 .048E-13 

5 .OOOr-07 

4 OOCE-OI 

2.S0E+0I 

1 .50E*01 

1392 

628 

0 

.00000 

1 .OOIE'OI 


.14400 

.11400 

1 .161E-I3 

2.055E-13 

1 .OOOE-Ol 

8. OOOE-Ol 

2.S0E^0l 

1 .50E*0l 

1155 

502 

0 

-00000 

1 .252E-01 


.11200 

.08711 

6.610E-14 

8 .593E-I4 

1 .250E-01 

1 .000F*00 

2.50E*0I 

1 .50E*OI 

987 

418 

0 

.00000 

1 .503E-01 


.11200 

.09073 

6.19IE-14 

1 .924E-13 

1 .500E-01 

1 . 200E*00 

2.50E*0I 

1.50E*0I 

809 

314 

0 

.00000 

2.001E 01 


.09700 

.07956 

6. 727E-1 3 

7.606E-13 

2.000E-01 

1 .600E*00 

2. 501*01 

1 .50E*0i 

683 

251 

0 

.00000 

2.503F-0I 


.08700 

.07228 

3.1 78Z-1 1 

3. 137E-I1 

2 .500E-01 

2.000E*00 

2 .50E*0I 

1 .50E*OI 

614 

209 

0 

-00000 

3.006E-01 


.07600 

.06277 

6.428E-10 

6.296E-10 

3.000E-01 

2.400E*00 

INTEGRATION METHOD: 

S5FSBD 






















MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLHSAL 

MINIMUM 

MAXIMUM 

ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP 

- TIME 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SUE 

2 .50E*01 

1 .t0E*01 

2590 

1256 

0 

.00000 

5.003E-02 


.22200 

.16618 

1 .576E-13 

1 .9#1E-13 

5 .OOOE-02 

3. OOOE-Ol 

2.50E*01 

1 -10E*01 

* 

628 

0 

.00000 

1 .OOlE-01 


.11200 

.08)01 

1 . 104E-13 

1 .655C-13 

1 .OOOE-Ol 

6. OOOE-Ol 

2.50E*0I 

1 .10E*0) 

1 104 

502 

0 

.00000 

1.2>2E-01 


.09650 

.07271 

1 .174E-13 

1 .272E-13 

1-250E-0I 

7.500E-01 

;.5or*oi 

1 .I0E*01 

936 

418 

0 

.00000 

1 .503E-0I 


.08000 

.05983 

1 .596E-13 

1 .986E-13 

1 .500E-01 

9. OOOE-Ol 

2.50E*01 

1 . 10E*0I 

739 

314 

0 

.00000 

2.001E-01 


.06800 

.05207 

7.206E-I2 

6 .915E-I2 

2.000E-01 

1 .200E*0D 

2 .50E*01 

1 .10E*0t 

624 

251 

0 

.00000 

2.503E-01 


.05200 

.03855 

1 .646E-10 

1 .611E-10 

2.500E-01 

1 .500C*00 

2.50E*01 

1 .10E*01 

551 

209 

0 

.00000 

3.006E-01 


.05000 

.0)812 

2.154E-09 

2. 14«E-09 

3. OOOE-Ol 

1 .800E*00 

INTEGRATION METHOD' 

ssFsri 






















MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ITER 

ORDER 

NFC 

NSTPA 

hSTPR 

NSTPT 

STEP SIZE 

CP 

- TIME 

OVNO 

ERROR 

ERROR 

STEP SIZE 

STEP SUE 

2.50C*Ot 

1 .10E*01 

2590 

1256 

0 

.00000 

5.003E-02 


.22300 

.16718 

I.576E-13 

1 .991E-13 

5.000E-02 

3. OOOE-Ol 

2.50C*01 

1 .I0E*01 

1345 

628 

0 

.00000 

1 .OOlE-01 


.11200 

.08301 

1 .104E-13 

1 .655E-13 

1 .OOOE-Ol 

6 . OOOE-Ol 

2.50E*0i 

1.10E*01 

1104 

502 

0 

.00000 

1 .252E-01 


.09450 

.07071 

1 .174E-13 

1 .272E-13 

1 .250C-01 

7.500E-0I 

2,50E*01 

1.I0E*01 

936 

418 

0 

.00000 

I.503E-01 


.08600 

.06583 

1 .596E-13 

1 .986E-i3 

1.500E-01 

9. OOOE-Ol 

2.50E*0l 

1 .10E*01 

719 

314 

0 

.00000 

2. 00 IE-01 


.06250 

.04657 

7.206E-12 

6.915E-12 

2. OOOE-Ol 

1 .200E*00 

2.50E*0I 

1.I0E«0I 

624 

251 

0 

.00000 

2.503E-01 


.05600 

.04255 

1 .646E-10 

1 .611E-10 

2.500E-01 

1 .500E*00 

2.501*01 

1 .10E*0I 

551 

20# 

0 

.00000 

3.0O6E-0I 


.05600 

.04412 

2.154E-09 

2.144E-09 

3. OOOE-Ol 

1 .800E*00 


Table IV. 1.1 

COMPAR Sunsnary of Statistics for the 
Harmonic Oscillator Problem 




r 


dt ■ *: 


!!i 

dt ■ *4 


djc, 

dt' 

"dT 


2 

3 


or by the set of second-order differential equations 


dx 

dt 


dx 

dt 


47 


where 


X.J *2 * 


For the initial conditions x^(0) = 1.0 , 


x_(0) = 0.0 , x-(0) = 0.0 , x.(0) = 1.0 , the analytic solution is 


given by: 


x^(t) = cos(t) 
X 2 <t) = f.in(t) 
x^Ct) = -sin(t) 
Xjj(t) = cos(t) 


The circular two-body orbit is shown in Figure IV. 1.2a along with an 
elliptic two-body orbit of eccentricity 0.6. 


By varying the final time of integration to 2 tt, 12ir and 20ir, 
the optimum orders for ABFS, KSGFS, SSFSBD and SSFSFE were found to be 
8, 13i Ilf respectively, for a range of stepsizes of 0.1 to 0.35. 
The efficiency curves for the variable-step integrators ODE, RK(7)8, 
RKN7(8), KRGH1 and KR0GH2 are shown in Figures IV. 1.2b through IV.1.2e. 
The efficiency curves for ABFS, KSGFS, SSFSBD and SSFSFE relative to 
those for ODE, KR0GH1 and KROGH2 are shown in Figures IV.1.2f through 






Figures IV. 1,2 b and c 

Efficiency Curves for the Circular Two Body Problem 





Efficiency Curves for the Circular Two Body Problem 






Figures IV. 1.2 f and g 

Efficiency Curves for the Circular Two Body Problem 
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Figures IV. 1.2 h and i 

Efficiency Curves for the Circular Two Body Problem 
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Figures IV. 1.2 j and k 

Efficiency Curves for the Circular Two Body Problem 




,.oi i-oi s-oi »-oi s-oi 9-01 4-ot 9-01 j-oioj-dJn-Oiji-dt 
U0UU3 maol') wnMixtjM 


S3 


ABSEKR 

RELERR 

NKF 

NS7PA 

NSrPR 

NSTPR/ 

NSTl 

AVERAGE 
lEP SIZE 

CP 

- TIME 

OVHD 

END POINT 
ERROR 

GLOBAL 

ERROR 

MINIMUM 
STEP SIZE 

MAX I MUM 
STEF MZE 

1 .OOE-OA 

I .OOE-12 

f»90 

53 

0 

.00000 

1 .186E*00 


.07200 

.03733 

2.352E-02 

2.3'2E-02 

2.000E-02 

1 .334E*00 

1 .00E«06 

1 -OOE-12 

1138 

89 

0 

.00000 

7.060E-01 


.11100 

.08641 

9.399E-04 

9.399E-04 

2.000E-02 

7.302E-01 

1 .OOE-08 

1 .OOE-12 

2016 

133 

0 

.00000 

4.054E-01 


.19100 

.14819 

2.169E-03 

2.169E-03 

2.000E-02 

4.261E-01 

1 .OOE-10 

1 .OOE-12 

3324 

271 

0 

.00000 

2.319E-01 


.32600 

.23116 

4.180E-07 

4. 180E-07 

2.000E-02 

2.406E-01 

l.OOE-n 

1 .OOE-13 

4694 

361 

0 

.00000 

1 740E-01 


.43800 

.33831 

3.668E-08 

3.668E-08 

2.000E-02 

1.803E-01 

i.ooE-n 

1. OOE-13 

6241 

480 

0 

.00000 

1 .309E-01 


.38300 

.43046 

7.646E-09 

7 .646E-09 

2 .OOOE-02 

1 .355E-01 

! .OOE-13 

1 .OOE-U 

8308 

639 

0 

.00000 

9.833E-02 


.77300 

.59656 

1 .031E-09 

1 .051E-09 

2.000E-02 

1.017E-01 

IHTEGRATION HETHOD. 

RK(7)8 






















MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ABSEU 

RELERR 

HFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP 

- TIME 

CVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

1 .OOE-04 

1 .OOE-12 

911 

70 

0 

.00000 

8.976E-01 


.05550 

.03613 

6.770E-02 

6.770E-02 

2.000E-02 

1 .OOIE^OO 

1 .OOE-06 

l.OOE-12 

1483 

lU 

0 

.00000 

5.512E-01 


.09600 

.06431 

8.121E-04 

8.121 E-04 

2.000E-02 

5.808E-01 

1 .OOE’Oa 

1 .OOE-12 

2349 

196 

0 

.00000 

3.206E-01 


.13800 

.10387 

5.697E-06 

3.697E-06 

2.000E-02 

3.338E-01 

1 .OOE-10 

l.OOE-12 

4486 

343 

0 

.00000 

1 .821C-01 


.26200 

.18673 

3.373E-08 

3.573E-08 

:,000E-02 

1 .884E-0I 

I.OOE-II 

1. OOE-13 

3969 

439 

0 

.00000 

1.369E-01 


.38700 

.26024 

2.861E-09 

2.861E-09 

2.000E-02 

1 .414E-01 

1 .OOE-I2 

1 .OOE-13 

7880 

606 

0 

.00000 

1 .037E-01 


.30300 

.33563 

2 .424E-10 

2 .424E-10 

2.000E-02 

1 .068E-01 

1 .OOE' :3 

1 .OOC-U 

9778 

732 

0 

.00000 

8.333E-02 


.62200 

.41434 

2.829E-U 

2.830E-!! 

2 .OOOE-02 

8.307E-02 

INTECRATIOII METHOD: 

ODE 























MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ABSERI 

RELEU 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP 

- TIME 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

1 .OOE-04 

l.OOE-12 

373 

186 

0 

.00000 

3.378E-01 


.09300 

.08308 

3.090E-01 

3.083E-OI 

2 .102E-03 

3.382E-01 

1 .OOE-06 

l.OOE-12 

609 

304 

0 

.00000 

2.067E-01 


.14600 

.13307 

2.939E-04 

3.167E-04 

2.102E-04 

2.133E-01 

1 .OOE<08 

l.OOE-12 

771 

385 

0 

.00000 

1 .632E-01 


.19300 

. 1 7663 

1 .472E-03 

4t"E-03 

2.102E-03 

1.722E-01 

I .OOE-IO 

l.OOE-12 

1071 

335 

0 

.000C3 

1 -174E-01 


.28400 

.26126 

1 .082E-08 

1 .098E-08 

2 .102E-06 

1 .378E-01 

1 .OOE-ll 

1 .OOE-13 

1301 

730 

0 

.00000 

8.378E-02 


.41200 

.38012 

4.223E-09 

9.402E-07 

6.648E-07 

8.7I3E-02 

1 .OOE-12 

1 .OOE-13 

1403 

701 

0 

.00000 

8.963E-02 


.38400 

.33420 

4.710E-09 

8.919E-07 

2 .102E-07 

1 .102E-01 

1 .OOE-13 

1. OOE-U 

2273 

1136 

0 

.00000 

3.331E-02 


.62000 

.37173 

3.419E-11 

1 .410E-06 

6 .648E-08 

6.971C-02 

UfTECllATlON METHOD: 

ER0CH2 






















MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ABSERR 

ULERR 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP 

- TIME 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

1 .OOE-04 

l.OOE-12 

397 

211 

0 

.00000 

2.978E-01 


.11100 

.10257 

2.001E-04 

1 .978C-04 

1 .OOOE-02 

3.200E-01 

I .OOE-06 

l.OOE-12 

426 

219 

0 

.00000 

2.869C-01 


.13700 

.12795 

1 .309E-03 

1 .293E-05 

1 .230E-03 

3.200E-01 

1 .OOE-08 

l.OOE-12 

810 

418 

0 

.00000 

1.303E-01 


.23400 

.23680 

8.084E-09 

8.074E-09 

7.813E-05 

1 .600E-01 

1 .00E<10 

l.OOE-12 

830 

434 

0 

.00000 

1 .448E-01 


.26800 

.23037 

1 .lOOE-IO 

1 .093E-10 

9.766E-06 

1 .600E-01 

l.OOE'll 

1. OOE-13 

831 

442 

0 

.00000 

1 .422E-01 


.27000 

.23193 

2.137E-10 

2.I45E-10 

2.441E-06 

1 .600E-01 

\ .OOE-12 

1. OOE-13 

1628 

827 

0 

.00000 

7.398E-02 


.31200 

.47743 

6.0D4C-11 

3.994E-11 

1 .221E-06 

8.000E-02 

1 .OOE-13 

1 .OOE-U 

1607 

833 

0 

.00000 

7.325E-02 


.30400 

.46987 

6.304C-11 

8 632C-07 

3.052E-07 

8.000E-02 

IIITECRATION HETHOD: 

RROCHl 






















MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ABSEU 

ULERR 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP 

- TIME 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SUE 

1 .OOE-04 

l.OOE-12 

426 

220 

0 

.00000 

2.836E-01 


.14100 

.13195 

3.784E-02 

3.736E-02 

1 .OOOE-02 

3.200E-01 

1 .OOE-06 

l.OOE-12 

812 

411 

0 

.00000 

1.329E-01 


.27100 

.23376 

7.462E-06 

7.424E-06 

1.230E-03 

1 .600E-01 

l.OOE-08 

l.OOE-12 

828 

423 

'O 

.00000 

1 .483E-01 


.28400 

.26642 

4.932E-06 

4.923E-06 

7.813E-03 

1 .600E-01 

l.OOS-10 

l.OOE-12 

1613 

818 

0 

.00000 

7.681E-02 


.35300 

.31874 

3.332E-10 

3.34IE-10 

9.766E-06 

8.000E-02 

1. OOE-ll 

1. OOE-13 

1623 

823 

0 

.00000 

7.616E-02 


.36600 

.33333 

4.223E-10 

4.210E-10 

2.441E-06 

8.000E-02 

l.OOE-12 

1. OOE-13 

1633 

829 

0 

.00000 

7.379E-02 


.36900 

.33432 

4.236E-10 

4.233E-10 

1 .221C-06 

8.000E-02 

1. OOE-13 

1. OOE-U 

2334 

1284 

0 

.00000 

4.893E-02 


.88400 

.83019 

2.036E-10 

8.632E-07 

3.032E-07 

8.000E-02 


Table IV. 1.2 

COMPAR Summary of Statistics for the 
Circular Two Body Problem 
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INTtCKATlON METHOD ABFS 


ITER 

ORDER 

NFE 

NSTPA 

NSTFR 

NSTPR/ 

NSTFT 

AVERAGE 
STEP SIZE 

CP - 1i ML 

2.50E^01 

8-OOE^OO 

1385 

628 

0 

.00000 

1 -OOlE-01 

.14500 

2.50E*01 

8 .00E*00 

1141 

502 

0 

.00000 

1 -252E-01 

.10700 

2.50E*01 

S.OOE^OO 

981 

418 

0 

.00000 

1 .503E-01 

.11300 

2.50Et01 

S.OOE^OO 

789 

314 

0 

OOOOO 

2.001E-01 

.07950 

2.50E+01 

8 .OOE^OO 

687 

251 

0 

.00000 

2.503E-01 

.07200 

2.50E*0l 

a.OOE^OO 

627 

209 

0 

. 00000 

3.006E-01 

,06350 

2.50E^01 

8 .00E*00 

209 

4 

0 

.00000 

0. 

0.00000 

INTEGRATION METHOD: 

KSCFS 










NSTPR/ 

AVERAGE 


ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SUE 

CP - IIME 

2.50E^01 

1.30E401 

1374 

628 

0 

.00000 

I .OOlE-01 

.16000 

2.50E^01 

1 .30E^01 

1135 

502 

0 

.00000 

1 -252E-01 

. 1 2800 

2.50E»01 

1 .30E^01 

967 

418 

0 

.00000 

1 .503E-01 

.11200 

2.50E^01 

1 .3OE»0! 

772 

314 

0 

.00000 

2 .OOlE-01 

.10400 

2.50E+01 

1.30E>0t 

672 

251 

0 

.00000 

2.503E-01 

. 10800 

2,50E^01 

1 .30E^01 

339 

6 

0 

.00000 

0. 

0.00000 

2 .50E^01 

1 .30Et01 

554 

179 

0 

.00000 

3.510E-01 

.08000 

INTEGRATION METHOD: 

SSFSBD 










NSTPR/ 

AVERAGE 


ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SUE 

CP - TIME 

2.50E^0l 

1 .lOE^Ol 

1356 

628 

0 

.00000 

1 .OOlE-01 

.14300 

2,50Et01 

l.lOE^Ol 

1104 

502 

0 

.00000 

1 .:52E-01 

.12700 

2.50E*01 

1 . tOE^Ol 

947 

418 

0 

.00000 

1 .503E-01 

.10900 

2.50E»0l 

1 .lOE^Ol 

750 

314 

0 

.00000 

2.001E-U1 

.08800 

2.50E*01 

1 .lOE^Ol 

646 

251 

0 

.00000 

2.503E-01 

.07700 

2.50Et0i 

l.lOE^Ol 

573 

209 

0 

.00000 

3.r06E-01 

.07150 

2.50E^0l 

1 . lOE^OI 

524 

179 

0 

.00000 

3.510E-01 

.06950 


MAXlMim 

END fOlNl M 'ML HINIML'M MAXIMUM 
OVHD ERKUK LLKvK SThT SUE STEP SUE 


.11559 i.I74E-06 I.W2E-06 1 .OOOE-Ol 5.000E-01 
.08277 8.208E-06 8.186E-06 1.250E-0! 6.250E-0I 
.09217 3.913E-05 3.897E-05 1.500E-0I 7.500E-01 
.06274 4.I97E-04 4 193E-04 2. OOOE-Ol I .OOOE^OO 
.05741 2.249E-03 2.243E-03 2.5O0E-OI 1 .250E*00 
.05018 1.039E-02 1.0I2E-02 3. OOOE-Ol i.500E*00 
0.00000 *** METHOD FAILED TO REACH FINAL TIME 


MAXIMUM 

END POINT GLOBAL MINIMUM MAXIMUM 
OVHD ERROR ERROR STEP SIZE STEP SIZE 


.13082 7.359E-13 2.222E-12 1 -OOOE-Ol 7. OOOE-Ol 
.10390 2.194E-11 2.I89E-I1 1.250E-0! 8.750F-01 
.09146 1.2245:-lO I.220E-10 I -500E-01 1 050E*00 
.08760 5.565E-09 5.56IE-09 2. OOOE-Ol 1 -400E^00 
.09373 9.680E-08 9.659E-08 2.500E-01 l.750E^00 
0.00000 METHOD FAILED TO REACH FINAL TIME 
.06823 1.972E-06 I -953E-06 3.500E-01 2.450E^00 


MAXIMUM 

END POINT GLOBAL MINIMUM MAXIMUM 
OVHD ERROR ERROR STEP SUE STEP SIZE 


.1 1420 3.078EM 1 3.076E-11 1 .OOOE-Ol 
.10355 2.674E-I1 2.670E-11 1..50E-01 
.08889 4.395E-11 4.389E-M I.500E-0I 
.07207 1.965E-09 1 .968E-09 2. OOOE-Ol 
.06328 4.388E-08 4.386E-08 2.500E-01 
.05933 4.869E-07 4.859E-07 3. OOOE-Ol 
.05837 3.028E-06 3.053E-06 3.500E-01 


6. OOOE-Ol 
7.5O0E-OI 
9. OOOE-Ol 
1 .200E^O0 
1.500E+00 
1 .SOOE^OO 
2.100E^O0 


INTEGRATION METHOn; SSFSFC 
ITER ORDER NFE NSTPA 

NSTPR 

NSTPR/ 

NSTPT 

AVERAGE 

STEP SUE CP - TIME 

OVHD 

END POINT 
ERROR 

MAXIMUM 

GLOBAL 

ERROR 

MINIMUM 
STEP SUE 

MAXIMUM 
STEP SUE 

2.50E^0I 

1 .10E«01 

1356 

628 

0 

.00000 

1 .OOlE-01 

.14300 

.11420 

3.078E-II 

3.076E-1I 

l.OOOE-OI 

6. OOOE-Ol 

2.50E«01 

I.10E«0t 

1104 

502 

0 

.00000 

1 .252E-01 

.12700 

.10355 

2.674E-11 

2.670E-I1 

1.250E-01 

7.500E-01 

2.50E*0I 

I.lOE^Ol 

947 

418 

0 

.00000 

1 .503E-01 

.11200 

.09189 

4.395E-U 

4.389E-11 

I.500E-01 

9. OOOE-Ol 

2.50E«01 

l.lOEtOl 

750 

314 

0 

.00000 

2.00IE-01 

.08850 

.07257 

1.965E-09 

I.968E-09 

2. OOOE-Ol 

1 . 20OE^OO 

2.50E«0I 

l.lOE^Ol 

646 

251 

0 

.00000 

2.503E-01 

.07950 

.06578 

4.388E-08 

4.386E-08 

2.5O0E-OI 

1.500E^OO 

2.50E^0I 

l.lOE^Ol 

573 

209 

0 

-OOOOO 

3.006E-01 

.07100 

.05883 

4.869E-07 

4.859E-07 

3.000B-0I 

l.BOOE^OO 

2.50E^0I 

l.lOE^Ol 

524 

179 

0 

.00000 

3.5I0E-01 

.06900 

.05787 

3.028E-06 

3.053E-06 

3.500E-01 

2.I00E«^00 


Table IV. 1,2 




COMPAR Summary of Statistics for the 
Circular Two Body Problem 
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IV.1,2i. Figures IV.1.2j through IV. 1.2m show the efficiency curves 
for ABFS, KSGFS, SSF5BD and SSFSFE relative to those of RK(7)8 and 
RKN7(8). To illustrate the efficiency curves, the integration interval 
was taken to be from 0.0 to 20 tt. A summary of these results are given 
in Table IV. 1.2, 


IV. 1 .3 Elliptic Problem of Two Bodies 


The elliptic problem of two bodies is modeled in two dimen- 
sions by the set of first-order differential equations 






3 


or by the set of second-order differential equations 


d^x. 




2 2 2 

where r s ♦ Xg . The analytic solution of this problem may be 
written in terms of the eccentricity, e , of the orbit as 


x^(t) s cos u 

XgCt) = (1-e^)^^^ sin u 

Xj(t) s -sin u/(1-e cos u) 

2 1/2 

Xjj(t) = (1-e ) cos u/(1-e cos u) 




wnere u is found by solving Kepler's equation 
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M = u - e sin u 

and M=t . The eccentricity was chosen to be equal to 0.6, and the 
initial conditions were 

x^(0)=0.4. X 2 ( 0 )= 0 . 0 , x^(0)=0.j and x^(0)=2.0 . 

Figure IV. 1.2a shows the eccentric two-body orbit relative to the cir- 
cular two-body orbit. 

In order to determine the optimum order for ABFS, KSGFS, 
SSFSBD and SSFSFE, the final time of the integration intervals examined 
were the same as those used in the circular two-body problem, i.e., 

12-rr and 20ti. However, due to the relatively large eccentricity of the 
orbit, the range of stepsize examined was from 0.001 to 0.025, about 
one-tenth of the values used for the circular two-body problem. The 
optimum orders were determined to be 12 for ABFS, 12 for KSGFS and 14 
for SSFSBD and SSFSFE. 

The efficiency curves for the variable-step integrators are 
shown in Figures IV. 1.3a through IV. 1.3d. The efficiency curves for 
ABFS. KSGFS, SSFSBD and SSFSFE relative to those for ODE, KR0GH1 and 
KR0GH2 are shown in Figures IV.1.3e through IV.1.3h. Figures IV.2.3i 
through IV. 2. 31 show the efficiency curves for ABFS, KSGFS, SSFSBD and 
SSFSFE relative to those of RK(7>8 and RKN7(8). The integration inter- 
val used to illustrate the efficiency curves was taken to be from 0.0 
to 20tt. The data for these figures are summarized in Table IV. 1.3. 




Figures IV. 1.3 a and b 

Efficiency Curves for the Eccentric Two Body Problem 
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Figures IV. 1,3 c and d 

Efficiency Curves for the Eccentric Two Body Problem 




Figures IV, 1,3 e and f 

Efficiency Curves for the Eccentric Two Body Problem 
















Efficiency Curves for the Eccentric Two Body Problem 
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i 


IKrtAV'ATION METH(JU 

N 

7(8) 











NSTPR/ 

AVERAGE 


END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

\B^LRR RELERR 

NFE 

NSTPA 

NSTPR NSTPT 

SrLF SIZE CP - TIME 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 


1 .OOE-OA 

1 .OOE-12 

117 

7 


.22^22 

6.967E-01 

0.00000 

0.00000 

*** METHOD FAILED TO REACH FINAL TIME 

1 .OOE-06 

1 .OOE-12 

2146 

121 

44 

.26667 

5.193E-01 

.20500 

.16247 

1.741E-01 

1.241E-01 

2.000E-02 

1.369E*00 

1 .OOE- 06 

1 .OOE-12 

3173 

109 

55 

.22541 

3.324E-01 

.29500 

.23212 

1 .670E-03 

1 .670E-03 

2.000E-02 

7.985E-01 

1 -OOE-IO 

1 .OOE-12 

4005 

30* 

7 

.02273 

2.087E-01 

.37500 

.29563 

1 .410E-05 

1 .410E-05 

2.000E-02 

4.515E-01 

1 .OOE-U 

1 .OOE-13 

5214 

401 

0 

.00000 

1 .567E-01 

.49000 

.38668 

1 .42IE-06 

1 .4 21E-06 

2.000E-02 

3.391E-01 

1 -OOE- 12 

1 .OOE-13 

6956 

535 

0 

.OoCOO 

1 .174E-01 

.66100 

.52316 

1 .197E-07 

I.197E-07 

1 .663E-02 

2.544E-01 

1 -OOE- 13 

1 .OOb-14 

9283 

714 

0 

-OOtfOO 

P.800E-02 

.86100 

.67704 

9.639E-09 

9.639E-09 

1 .231E-02 

1 .909E-01 

INTEGRATION METHOD: 

RX(7: 

)8 




















MAXIMUM 








NSTPR/ 

AVERAGE 



END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ABSERR 

RELERR 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP - TIKE 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

1 .OOE-04 

1 .OOE-12 

1792 

111 

29 

.20714 

5.661E-01 

.10000 

.06449 

2.087E-0I 

2.087E-01 

2.000E-02 

l.771E*00 

1 .OOE-06 

1 .OOE-12 

2852 

163 

61 

.27232 

3.855E-01 

.17700 

.12048 

4.506E-03 

4.506E-03 

2.000E-02 

1 .097E+00 

1 .OOE-08 

1 .OOE-12 

4597 

264 

97 

.26870 

2.380E-01 

.28400 

.19290 

3.361E-05 

3.361E-05 

2.000E-02 

6.642E-01 

1 OOE-IO 

1 .OOE-12 

757C 

450 

144 

,24242 

1.396E-01 

.46901 

.31881 

1.864E-07 

1.864E-07 

2.000E-02 

3.783E-01 

\ .OOE-1 1 

1. OOE-13 

8854 

584 

105 

.15239 

1.076E-01 

.56300 

.38754 

7.766E-09 

7.766E-09 

2 .OOOE-02 

2.843E-01 

1 .OOE-12 

1 .OOE-13 

9778 

752 

0 

.00000 

8.355E-02 

.61600 

.42723 

3.751E-09 

3.751E-09 

1 .681E-02 

2.171E-01 

1 .OOE- 13 

1 .OOE-14 

12000 

922 

1 

.00108 

6.815E-02 

.75800 

.52020 4.842E-10 

4.842E-10 

1.2lOE-(u 

1 .805E-01 

INTEGRATION METHOD: 

ODE 





















MAXIMUM 








NSTPR/ 

AVERAGE 



END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ABSERR 

RELERR 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP - TIME 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

I .OOE-04 

1 .OOE-12 

1152 

558 

0 

.00000 

1.126E-01 

.32100 

,29817 

1.736E+00 

2.508E*00 

9.759E-04 

4.997E-01 

1 .OOE-06 

1 .OOE-12 

1776 

870 

0 

,00000 

7.222E-02 

.51400 

.47881 

1 .996E-02 

1 -996E-02 

9.759E-05 

2.599E-01 

1 .OOE-08 

1. OOE-12 

2644 

1304 

0 

.00000 

4.818E-02 

,78700 

,73461 

1 .820E-04 

1.819E-04 

9.759E-06 

1.837E-01 

1 .OOE-IO 

1 .OOE-12 

3666 

1819 

0 

.00000 

3.454F-02 

1 .10100 

1 ,02835 

6.909E-07 

6.404E-06 

9.759E-07 

9.533E-02 

1 .OOE-n 

1 .OOE-13 

4374 

2175 

0 

.00000 

2.889E-02 

1.31100 

1 ,22432 

3.779E-08 

6.076E-CO 

3.086E-07 

8.2?;e-02 

I .COE- 12 

1 .OOE-13 

5462 

2720 

0 

.00000 

2.310E-02 

1 .64700 

1,53876 

1.172E-09 

4.483E-06 

9.759E-08 

7.460E-02 

1 .OOE- 13 

1 .OOE-14 

7055 

3522 

0 

.00000 

1 .784E-02 

2.05800 

1.91819 

2.448E-09 

6.278E-06 

3.086E-08 

6.671E-02 

INTEGRATION METHOD: 

KR0CH2 




















MAXIMUM 








NSTPR/ 

AVERAGE 



END POINT 

GLOBAL 

MINIMUM 

KAXlHim 

ABSERR 

:lerr 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP - TIME 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

1 .OOE-04 

1 .O^L-12 

1699 

844 

0 

.00000 

7.445E-02 

.51200 

.47833 

1 .041E-01 

1 .041E-01 

2.500E-03 

3.200E-01 

1 .OOE-06 

1 .OOE-12 

2429 

1214 

0 

.00000 

5.I76E-02 

.73800 

.60987 

1.938E-04 

1 .917E-04 

1.563E-04 

1 .600E-01 

1 .OOE-08 

1 .OOE-12 

3461 

1723 

0 

.00000 

3.647E-02 

1 .07100 

1.00242 

1.179E-05 

1.179E-05 

1 .953E-05 

1.600E-01 

1 .OOE-IO 

l.OOE-l? 

4400 

2197 

0 

.00000 

2.860E-02 

1.38300 

1.29581 

5.119E-08 

5.109E-08 

2.441E-06 

8.000C-02 

1 .OOE-U 

1 .OOE-13 

4976 

2496 

0 

.00000 

2.517E-02 

1 .56400 

1.4 39 

1.3l9F>)8 

4.005E-06 

6.I04E-07 

8.000E-02 

1 .OOE-12 

1 .OOE-13 

6339 

3162 

0 

.00000 

1 .987E-02 

1 .99800 

1.8/238 

5.318E-10 

4.005E-06 

1 .526E-07 

8.000E-02 

1 .OOE- 13 

1 .OOE-14 

7442 

3712 

0 

,00000 

1.693E-02 

2.32300 

2.17553 

1.489E-10 

4.005E-06 

7.629E-O0 

8.000E-02 

INTEGRATION METHOD: 

KROGHl 




















HAXINim 








NSTPR/ 

AVERAGE 



END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ABSERR 

RELERR 

NFE 

HI TP^ 

NSTPR 

NSTPT 

STEP SIZE 

CP - TIKE 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

1 .OOE-04 

1 .OOE-12 

1729 

849 

0 

.00000 

7.401E-02 

.58700 

.55274 

2.125E-01 

2.133E-01 

1.250E-03 

3.200E-01 

1 .OOE-06 

1. OOE-12 

2448 

1213 

0 

.00000 

5.180E-02 

.83600 

.78749 

1 -489E-03 

1.477E-03 

1.563E-04 

1.600E-01 

1 .OOE-08 

1. OOE-12 

3624 

1810 

0 

.00000 

3.471E-02 

1 .26500 

1.19318 

1.624E-05 

1.624E-05 

I.953E-05 

l,600E-0i 

1 .OOE-IO 

1. OOE-12 

4427 

2206 

0 

.00000 

2.848C i2 

1.57100 

1.48327 

I.U4E-07 

I.I12E-07 

2.441E-06 

8.000E-02 

1 .OOE-11 

1. OOE-13 

5104 

2544 

0 

.00000 

2.47wE-0: 

1.83100 

1.72986 

3.100E-09 

4.005E-06 

6.104E-07 

8.000C-02 

1 .OOE-12 

1. OOE-13 

6588 

3286 

0 

.00000 

1.912C-02 

2.35900 

2.22845 4.053E-10 

4.005E-06 

1.526C-07 

8.000E-02 

1 .OOE- 13 

1. OOE-14 

8249 

4124 

0 

.00000 

I.524B-02 

2.89500 

2.7315i 

1.722E-10 4.005E-06 

7.629C-08 

4.000E-02 


* Table IV. 1.3 

COMPAR Summary of Statistics for the 
Eccentric Two Body Problem 


“'V 
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INIm.kATION TMOn ARFS 


N'.ri'R/ AVtKACF 

ITER ORrEP nee NS TEA NSTER N'^TPT STEP SU.E CP 



1 .20K*0I 

1689; 

8376 

0 

.00000 

7.^0lE-03 1 

2 -MIMOI 

1 .20E*0l 

1 2 721 

6282 

0 

.00000 

1 .000E-U2 1 

7 .sor»(H 

1 .20E*0I 

10219 

S02S 

0 

. 01)000 

1.2SOE-02 1 

2 

1 .20t»0l 

8^43 

4187 

0 

.00000 

1 .SOlE-02 

2 .SUE* 01 

1 .20E»01 

6461 

3140 

0 

. 00000 

2.001E-02 

: soMOl 

1 20E*01 

S217 

2M2 

0 

.00000 

2.301E-02 

INTK.RATION HKTHOO 

users 










HSTPR/ 

AVERAGE 

UER 

ORDER 

NFE 

NSrPA 

NSTPR 

KSTPT 

STEP SUE CP 

2 .SOE*Ol 

1 .90E*0l 

16886 

8376 

0 

.00000 

7.SOU-03 1 

2 .Mir*oi 

\ .90E*0I 

12717 

6282 

0 

.00000 

l.OOOE-02 1 

2 .soFun 

1 .'»or.*oi 

10203 

S02S 

0 

.00000 

1.2S0E-02 ! 

2.SOF*Ol 

1 .90E*0i 

8S46 

4187 

0 

.00000 

l.SOlE-02 1 

2 .S0F*01 

I .«)0E*0! 

6471 

3140 

0 

.00000 

2.001E-02 

2 .SUF*Ol 

1 

S21S 

2M2 

0 

.00000 

2.S01E-02 


- TIME 

OVHD 

FND POINT 
ERROR 

MAXIMUM 

GIOBAL 

ERROR 

MINIMUM 
STEP SUE 

MAXIMUM 
STEP SUE 

1 .9^800 

1 .62 316 

3.110E<I1 

4.092E-11 

7 .300E-03 

3.2S0E-02 

1 .48^00 

1 .23291 

3.64U-10 

3 632E-10 

I .OOOE-02 

7 .OOOE'02 

1.16700 

.96430 

1 .326E-08 

1 . 32 7E-08 

1 230E-02 

8 . 7SOE-02 

.98000 

.81071 

1 .243E-07 

1 .241E-07 

1 .300E-02 

1 .030E-01 

.74700 

.61897 

4.046E-06 

4 .03SF.-06 

2 OOOE-02 

1 .400E-01 

.60900 

.S0362 

6.S66E-0S 

6 .336F-0S 

2 .SOOE-02 

1 .7S0E-01 





MAX 1 MUM 





END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

- TIME 

OVHD 

ERROR 

ERROR 

STEP SUE 

STEP SUE 

1 .98700 

1 .63238 

9.633E-U 

9.619E-11 

7 .300G-03 

7.300E>02 

1 .31300 

1 .26099 

7 .639E-10 

7 .633E~10 

I .OOOF. 02 

1 .OOOE-OI 

1 21100 

1 .00881 

7 .37 7E-10 

7.364E'lO 

1 .230E-02 

1 .230E-01 

1 .01800 

.84863 

1 .488E-09 

1 .4P2E 09 

1 .300E-02 

1 .300E-01 

.79800 

.66977 

6.46 3E-08 

6.471F 08 

2.000E-02 

2 .OOOE-Ol 

6 3300 

.33166 

3.898E-07 

3.976E-07 

2 .300E-02 

2 .SOOF-01 


INTH.RATION Mt THOD SS^SBD 

KAXIMUH 


ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPRy 

NSTPT 

AVERAGE 
STEP SUE CP 

> TIME 

OVHD 

END POINT 
ERROR 

GLOBAL 

ERROR 

MINIMUM 
STEP SUE 

MAXIMUM 
STEP SUE 

2.30E*01 

1 ,40E*0l 

16831 

8376 

0 

. 00000 

7.301E-03 I 

. 79600 

1.46207 

4.176E-11 

4.307E-11 

7 .300E-03 

6.000E-02 

2 .S0FM)1 

1 .40E*01 

12677 

6282 

0 

.00000 

l,O00E-02 I 

1-37200 

1 .121. '9 

6.413E-11 

6.464E-1 1 

1 .OOOC-02 

B.OOOE'02 

2 .30E*0l 

1 .40F*0l 

10163 

3023 

0 

.00000 

1 .230E-02 1 

.09200 

.89060 

4.131E-11 

4.293E-1 1 

1 .230E-02 

1 .OOOE-Ol 

2 ,3UF,*01 

1 .40E*01 

8487 

4187 

0 

.00000 

1 .S0ir.-02 

,90400 

.73382 

1 .204E-10 

1 . 303E-10 

1 .300E>02 

1 .200E-01 

2.30F*01 

1 .40F*0l 

6407 

3140 

0 

.00000 

2.001E-02 

.69600 

.36904 

2.827E-09 

4.782E-09 

2.000E-02 

1 .600B-01 

2 . 301*01 

1 .40E*01 

3131 

2312 

0 

.00000 

2.301E-02 

.34900 

.44693 

4.384F.-08 

1 .023E-07 

2 .SOOE-02 

2. OOOE-Ol 


INTEC.RATION HE I HOD SSKSPE 


N3TPR/ AVERAGE 


END POINT 


NAXlHim 

GLOBAL 


NINIMUN NAXlHim 


ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SUE 

CP - TIME 

OVHD 

ERROR 

ERROR 

STEP SUE 

STEP SUE 

2 .30E*01 

1 .40E*0l 

KS831 

8376 

0 

.00000 

7.301E-03 

1 .79100 

1 .43707 

4.176E-U 

4.307E-M 

.300E-03 

6.000E-02 

2 .30P*0l 

1 .40E*0l 

1 2677 

6282 

0 

.00000 

1 .OOOE-02 

1.33900 

1 .107 79 

6.413E-I 1 

6.464E-11 

1 .OOOE-02 

8.000E-02 

2.30E*01 

1 .40E*Ot 

10163 

3023 

0 

.00000 

1 .2S0E-02 

1 .09)00 

.89160 

4.131E-11 

4.293E-1 1 

1 ,230r.-02 

1 .OOOE-Ol 

2 .30E*01 

1 .40E*01 

8487 

4187 

0 

. 00000 

1 .301E-02 

.91700 

.74882 

1 .204E-10 

1 .303E-10 

1 .300E-02 

1 .700E-01 

2 .30E*0l 

1 .40E*01 

6407 

3140 

0 

. 00000 

2.001E-02 

.69900 

.37204 

2 .827E-09 

4.782E-09 

2 .OOOE-02 

1 .600E-01 

2.30E*0l 

1 .40E*0l 

3131 

2312 

0 

.00000 

2.301E-02 

.33900 

.43693 

4.384C-08 

1 .023E-07 

2 .300E-02 

2. OOOE-Ol 


Table IV. 1.3 

COMPAR Summary of Statistics for the 
Eccentric Two Body Problem 



IV. 1.^ Ful«'r Ripid-Body Problou 


The Euler equations of motion for a rit’.id body without exter- 
nal forces may be written as 


dx 


dt ■ *2*3 ‘ 


dx^ 

dt ' "*1*3 * 


dx^ 

'dt' ' • 


dx^ dx^ dx^ 

By defining Xjj = , x^ - -j- and x^ * jt ’ original set of 
three first-order ordinary differential equations may be extended to 
six equations. Tlie new set of six differential equations may be writ- 
ten in two forms. The first form is given by the set of first-order 
differential equations 


dx 


dt 

dX2 

'dt“ 

dx^ 

"dT 


*2*3 


■*1*3 


-.51x^X2 


dx^ 

"dT 


dt 


^*6 

"dt" 


-x^(x^ ♦ .51x2) 


-X2(x^ - .51x^) 


-.51x^(X2-x^) 


or the equivalent set of second-order differential equations 


d^x. 


dt 


2 2 
-x^(x^ ♦ *51 x2> 


d X, 


dt 


2 2 

- .51x^> 




dt 


2 2 
-.51x^(x 2 - x^ ) 



The second form is given by the set of first-order differential equa- 


dx 


dt ^ 
dx^ 


dt 

dt 


= X, 


= X, 


dx^ 

~dt 

dt 

dx^ 

dt 


X5X3 + *6*2 


- ( XjjX^ + X^ Xg ) 


-,51 (Xj^X^ + x^x^ ) 


or the equivalent set of second-order differential equations 


d^x, 


dt 


* 5*3 * 6*2 


d^x. 


dt 


-(x^x^ + *1*5^ 


d^x. 


dt 


-.51 (x^x^ x^x^ ) 


By expanding the state from three to six elements, the right-hand side 
of the differential equations may be written as a function of x^, 
and Xj or of x^, Xg. x^, x^, x^ or Xg. The -lytic solutions are 
given by the Jacobian elliptic functions 


x^ (t) = sn(t! . 51 ) 

X 2 (t) = cn(ti.51) 
x^(t) = dn(ti.5D 

with the initial conditions x^ 


x^(t) = cn(t! .51)dn(t! .51) 
xAt) s -sn(ti .51)dn(tl .51) 

D 

Xg(t) = -.51 sn(t| . 51 )cn(tl .51 ) 

s 0.0, x^(0) = 1.0, x^CO) = 1.0, 



X, (0) = 1.0, x,.(0) r 0.0 and x,(0) = 0.0. The Jacobian elliptic func- 

s D 

tions are poriolic. The functions sn(t|.51) and cn(t|.5D have a 
quarter-period of K and dn(tl.51) has a half-period of K where 
K = 1.8626 'i 080233273855203... The functions sn(t| .51 ) , cn(t| .51) and 
dn(t|.5D are shown in Figure IV. 1.4a. The integration interval was 
taken to be from 0.0 to 8K. 

To obtain the efficiency curves for the fixed-step/fixed- 
order integrators, the stepsize was varied from 0.03 to 0.10. The 
optimum orders for ABFS, KSGFS, SSFSBD and SSFSFE were 11, 16, 14 and 
14 for the first form of the differential equations and 11, 11, 10 and 
10 for the second form of the differential equations. 

The efficiency curves presented in Figures IV. 1.4b through y 
are for an integration interval from 0,0 to 8K. The efficiency curves 
for solving the first form of the differential equations are shown in 
Figures IV. 1.4b through IV. 1.4m. Figures IV. 1.4b through IV.1.4e show 
the efficiency curves for the variable-step integrators. Figures 
IV. 1.4f through IV.1.4i show the efficiency curves of ABFS, KSGFS, 
SSFSBD and SSFSFE relative to those for ODE, KR0GH1 and KR0GH2. Fig- 
ures IV.1.4J through IV, 1.4m show the efficiency curves of ABFS, KSGFS, 
SSFSBD and SSFSFE relative to those for RK7(8) and RKN7(8). Similar 
comparisons for the solution of the second form of the differential 
equations are shown in Figures IV.1.4n through IV. 1.4y. The results of 
the comparisons are summarized in Tables IV.1.4A and IV. 1.4B. 
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Figure IV. 1.4a 

Jacobian Elliptic Functions: 
sn(tl.51), cn(t|.51), and dn(tl.51) 
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Figures IV, 1.4 b and c 

Efficiency Curves for the Euler Rigid Body Problem 

F^F(T,Y) 





Efficiency Curves for the Euler Rigid Body Problem. 

F-F(T,Y) 





Figures IV. 1.4 f and g 

Efficiency Curves for the Euler Rigid Body Problem, 

F»F(T,Y) 






Efficiency Curves for the Euler Rigid Body Problem, 

P-F(T,Y) 





Figures IV. 1,4 j and k 

Efficiency Curves for the Euler Rigid Body Problem 

F-F(T,Y) 




Figures IV. 1. A n and o 

Efficiency Curves for the Euler Rigid Body Problem, 






M I71I 



Efftclcacy Curves for the Euler Rigid Sudy Problem 

F-F(T,Y,Y^^^) 
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Figures IV, 1.4 r and s 

Efficiency Curves for the Euler Rigid Body Problem, 

P"F(T,Y 
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Figures IV, 1,4 t and u 

Efficiency Curves for the Euler Rigid Body Problem 

F»«F(T,Y,Y^^^) 
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Figures IV, 1.4 v and w 

Efficiency Curves for the Euler Rigid Body Problem, 













S3 







N‘;ipR 

AVFRACK 




END POINT 

U LORAL 

MINIMUM 

MAXI Him 

KR 

KM tRP 

mf: 

NSU’A 

NSU'K 

sM 1 : 

<]rV SUF, 

cr 

- TIME 

OVHD 

ERROR 

ERROR 

.STEP SIZE 

STEP SUE 

1 .OOf-04 

1 .UOE'12 

196 

14 

1 

.06667 

1 064E*00 


.02100 

.01674 

5.173E-01 

5.173E-01 

2.000E-02 

1 .548E*00 

I OOE-06 

1 .OOE-12 

287 

22 

0 

.00000 

6.773E-01 


.03250 

.02626 

3.419E-03 

3.419E-03 

2.000E-02 

9.12IE-01 

1 .OOF-08 

1 .out -12 

482 

37 

0 

.00000 

4 .O27E-01 


.05600 

.0^552 

4 .804E-06 

4 .804E-06 

2 -OOOE-02 

5.337EOI 

1 OOE-10 

1 .OOE-12 

846 

65 

0 

.00000 

2.292E-OI 


.08800 

.06960 

5 - 392E-08 

5.392E-08 

2.000E-02 

3.019E01 

! -oin-i 1 

1 .oot-n 

1093 

84 

0 

.00000 

1 . 7 74E-U1 


.11100 

.08723 

5.120E-09 

5.120E-09 

2.000E-02 

2.261E-01 

\ .OOE-12 

1 .OOF.- 13 

1457 

112 

0 

.00000 

1 .330E-01 


.15800 

.12631 

4.779E-10 

4.779E-10 

2.000E-02 

1 .690E-0I 

1 -OOE-13 

1 .OOF-U 

1925 

148 

0 

.00000 

1 .007E-01 


.21700 

.17514 

4. 160E-1 1 

4.160E-11 

2 .OOOE-02 

1 .272E-01 

1 NTH. RATI UN MM HOD 

KK( 7)8 






















HAXIMUH 








NSTPR/ 

AVERACC 




END POINT 

C LORAL 

MINIMUM 

MAXlMim 

AA.SFRR 

REl.tRR 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SiZE 

CP 

- TIHE 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

1 .OOE'04 

1 .OOE-12 

232 

15 

3 

.16667 

9.93..E-01 


.01833 

.01329 

2.993E-03 

2.993E-03 

2 .OOOE-02 

1 .6 34E*00 

1 .OOE-06 

1 .OOE-12 

360 

2) 

5 

.17857 

6.479C-01 


.02550 

.01767 

6.948E-06 

6.948E-06 

2.000E-02 

9.304E-01 

\ .oor-08 

1 .OOE-12 

553 

36 

7 

.16279 

4. 139F-0I 


0406? 

.02864 

5 .994E-08 

5.994E-08 

2.000E-02 

6.003E-0I 

1 -OOE-iO 

1 .OOE-12 

852 

59 

7 

.10606 

2.526E-01 


.06600 

.04 74 7 

7.220E-10 

7.220E-10 

2.000E-02 

3.749E-01 

1 .OOE-ll 

1 .OOE-13 

1135 

77 

It 

.12500 

1 .935E-01 


.08750 

.06282 

8.914C-11 

6.914E-11 

2.000E-02 

3.026E-01 

1 .oor -12 

1 OOE-13 

U21 

99 

11 

.10000 

1 .505E-01 


.10500 

.07410 

1 .285E-1I 

1 .285E-11 

2.000E-02 

2.395E-01 

1 .OOK-13 

1 .OOE-14 

1 747 

125 

10 

.07407 

1 .192E-01 


. 1 3800 

.10001 

2.146E-12 

2.146E-12 

2.000E-02 

1 .981E-01 

INTK.RATION METHOD 

ODE 























HAXIMUH 








NSTPR/ 

AVERAGE 




END POINT 

GLORAL 

MINIMUM 

MAXIMUM 

ARSLRR 

RE I. ERR 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP 

- TIHE 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SUE 

1 .00E'04 

1 .OOE-12 

173 

86 

0 

.00000 

1 . 7 3 3E-01 


.05400 

.05024 

1 .170E-04 

1 .185E-04 

2 .039E-03 

2.610E-01 

1 .OUE'06 

1 .OOE-12 

111 

155 

0 

.00000 

9.6I4E-02 


.09400 

.05724 

4.464E-07 

4.865E-07 

2.039E-04 

1 .044E-01 

1 .OOF-08 

1 .OOE-12 

401 

700 

0 

.00000 

7 .4ME-02 


.12900 

.12028 

6. 723E-09 

7.408E-09 

2 .039E-05 

8.352E-02 

1 .OOE'IO 

1 OOE-12 

559 

2 79 

0 

.00000 

5.341E-02 


. 1 7600 

. 16384 

3.269E-11 

4.00IE-1 1 

2.039E-06 

6.681E-02 

\ .OOE'U 

1 OOE-13 

765 

382 

0 

.00000 

3.901E-02 


.25700 

.24036 

4.309E-12 

9.693E-07 

6.448E-07 

4.226E-02 

1 OOF' 12 

1 .OOE-13 

791 

395 

0 

.00000 

3. 7;:e-o2 


.27400 

25680 

5.I98E-I2 

9.196E-07 

2 .039E-07 

5.345E-02 

1 .00E'l3 

1 .OOE-14 

1055 

527 

0 

.00000 

2 .828E-02 


.35000 

. ’2706 

6.548E-12 

1 .4S4E-06 

6.448E-08 

3.380E-02 


INFFURAUUN METHOD 
ABSERR KELFRR 

krcx;h2 

NFE NSTPA 

NSTPR 

NSTPR/ 

NSTPT 

AVERAGE 
STEP SIZE 

CP - TIME 

OVHD 

END POINT 
ERROR 

MAXIMUM 

GLOBAL 

ERROR 

MINIMUM 
STEP SIZE 

MAXIMUM 
STEP SIZE 

1 .OOE-04 

1 .OOE-12 

245 

134 

0 

.00000 

1 .112E-01 

.09500 

.0896, 

2.681E-03 

2.681E-03 

1 .OOOE-02 

1 .600E-01 

1 .OOE'06 

1 .OOE-12 

382 

196 

0 

.00000 

7.603E-02 

.14500 

. 1 3669 

9.412E-06 

9.276E-06 

6.250E-04 

1 .600E-01 

1 .OOE-08 

1 .OOE-12 

523 

272 

0 

.00000 

5.478E-02 

. 20800 

. *9663 

2.117E-07 

2.117E-07 

7 .81 3E-05 

8.000E-02 

1 .OOE-IO 

1 .OOE-12 

705 

368 

0 

. 00000 

4.049E-02 

.77200 

,25667 

7.426E-09 

7.386E-09 

9. 766E-06 

6.000E-02 

1 OOE-M 

1 .OOE-13 

773 

409 

0 

.00000 

3.643B-02 

. 30200 

.28519 

1 .437E-12 

1 .429E-12 

2.441E-06 

4.000E-02 

1 .OOE-12 

1 .OOE-13 

818 

422 

0 

.00000 

3.531E-02 

.31900 

.30121 

7.159E-11 

9.176E-07 

6. 104E-07 

4.QOOE-02 

1 .OOE-13 

1 .UOE'14 

ISOT 

771 

0 

.00000 

1 .933E-02 

.57800 

.54518 

1 .176E-11 

9.176E-07 

3.052E-07 

4.000E-02 

INTEGRATION METHOD. 
ABSERR REI.ERR 

KROGHl 
NFE NSTPA 

NSTPR 

NSTPR/ 

NSTPT 

AVERAGE 
STEP SIZE 

CP - TIME 

OVHD 

END POINT 

ERROR 

MAXIMUM 

GLOBAL 

ERROH 

MINIMUM 
STEP SHE 

MAXIMUM 
STEP SIZE 

1 .OOE-04 

1. OOE-12 

284 

144 

0 

.00000 

1 .015E-01 

.12400 

.11782 

2.548E-05 

3.479E-05 

1 .OOOE-02 

1 .600E-01 

1 .OOE-06 

1 .OOB-12 

384 

201 

0 

.00000 

7.413E-02 

.17500 

.16665 

3.509E-07 

4 .069E-07 

6.'50E-04 

1 .600E-01 

1 .OOE-08 

1 .OOE-12 

531 

269 

0 

.00000 

5.539E-02 

.24200 

.23045 

1 .530E-09 

1 .555B-09 

7 3E-05 

8.000E-02 

1 .OOE-10 

1. OOB-12 

732 

371 

0 

.00000 

4.016E-02 

.33500 

,31908 

7.374C 10 

7.367E-10 

9./66E-06 

8.000E-02 

1 .OOE-ll 

1. OOE-13 

794 

406 

0 

.00000 

3.670E-02 

.37500 

.35773 

4.579K-12 

4.552E-12 

2.441E-06 

4.000E-02 

1. OOE-12 

1. OOE-13 

1406 

714 

0 

.00000 

2.087B-02 

.65500 

.62442 

3.470E-12 

9.I76E-07 

6.104E-07 

4 .OOOE-02 

1 .OOE-13 

1 .OOE-14 

1542 

783 

0 

.00000 

1 .903E-02 

.70200 

.66847 

7.182E-12 

9.176E-07 

3.0'2B-07 

2 .OOOE-02 


Table IV.l.AA 


COMPAR Summary of Statistics for tbe 
Eulei Rigid Body Problem F*F(T,Y) 



lNrH.RATU>N METHOD- ABFS 


MAXIMUM 







NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

HAXIMI^ 

ITER 

ORDER 

KFE 

NSTPA 

NSTPR 

HSTPT 

STEP SUE 

CP 

- TIKE 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

2.'iOE*OI 

I .ICE*01 

1 lU 

496 

0 

.00000 

3.004E-02 


.16800 

.14377 

1 .670E-13 

2 560E-13 

3.000E-02 

1 .800E-01 

2.50E*01 

1 .lOE^Ol 

888 

372 

0 

.00000 

it.006E-02 


.12800 

.10869 

6.615E-13 

7.190E-13 

4.000E-02 

2.400E-01 

2.S0E+01 

1 .lOE+01 

740 

298 

0 

.00000 

5.000E-02 


.12500 

.10891 

5.530E 12 

6.121E‘12 

5 -OOOE-02 

3.000E-01 

2.S0E*01 

1 .lOE^Ol 

640 

248 

0 

.00000 

6.009E-02 


.09600 

.08208 

5.403C-11 

5.527E-11 

6.000E-02 

3.600E-01 

2.S0E*01 

1 .lOE^Ol 

579 

212 

0 

.00000 

7.029E-02 


.09550 

.08291 

3.620E-10 

3.234E-10 

7.000E-02 

4.200E-01 

2.50E*01 

1 .10E*01 

527 

186 

0 

.00000 

8.011E-02 


.07750 

.06604 

1 .385E-09 

1.471E-09 

8.000E-02 

4.800E-01 

2.S0E»01 

1 .lOE^Ol 

496 

165 

0 

.00000 

9.031E-02 


.07900 

.06821 

5.724E-09 

5.421E-09 

9.000E-02 

5.400E-01 

2.50E*01 

1 .I0E*0I 

464 

149 

0 

.00000 

1 .OOOE-01 


.07900 

.06891 

3.165E-08 

3.090E-0C 

1 -OOOE-01 

6.000E-01 

INTEGRATION METHOD; 

RSCFS 






















MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MlNlHim 

MAXIMUM 

ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP 

- TIME 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

2.SOE^Ol 

I .60E+01 

IM9 

495 

0 

.00000 

3.010E’02 


.19600 

.17166 

1.935E-11 

1.929E-11 

3.000C-02 

2.700E-01 

2.M)E*0l 

1 .60E*01 

871 

371 

0 

.00000 

4.016E-02 


.15600 

.13.06 

1.197E-11 

1.199E-11 

4.000E-02 

3.600E*01 

2 .^OE*Ol 

1 .60E*0l 

739 

297 

0 

.00000 

5.oi;e-o2 


.14400 

.12793 

4.135E-ll 

4.128E-11 

5 .OOOE-02 

4.500E-01 

2.^0E*0I 

1 .60E»0I 

639 

247 

0 

.00000 

6.033E-02 


.12500 

.11110 

2.180E-II 

2.176E-11 

6.000E-02 

5.400E-01 

2.S0E*0I 

I .60E+01 

583 

211 

0 

.00000 

7.062E-0? 


. 1 2800 

.11532 

3.032E-10 

2.990E-10 

7.000E-02 

6.300E-01 

2.50E*01 

1 .60E*^0l 

515 

185 

0 

.00000 

8.055E~02 


.11200 

.10080 

2 -666E-09 

:.638E-09 

8.000E-02 

7.200E-01 

2.^0E*01 

1 .60E+01 

505 

164 

0 

.00000 

9.086E-02 


.11200 

.10102 

1 .491E>08 

1 .468E-08 

9.000E-02 

8.100C-01 

2.^0E401 

1 -60b*0! 

457 

148 

0 

.00000 

1 -00 7E-01 


.11500 

,10506 

5.953E-08 

5.880E-08 

1 .OOOE-01 

9.000E-01 

INTEGRATION i 'THOD: 

SSFSBD 






















MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP 

- TIME 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 


2.50E»01 

1 .40E«01 

1089 

495 

0 

.00000 

3.010E-02 

.17300 

.14932 

5.526E-12 

5.529E-12 

3.000E-02 

2.400E-01 

2.50E*0l 

1 .40E*01 

85S 

371 

0 

.00000 

4.016E-02 

.14300 

.12441 

6.780E-12 

6.789C-12 

4 .OOOE-02 

3.200E-01 

2 .5OE*01 

1 .40E^01 

721 

297 

0 

.00000 

5 .01 7E-02 

,11300 

.09732 

7.145E-12 

7.221E-12 

5 -OOOE-02 

4 .OOOE-Ot 

2.50E^01 

1 .40E«01 

621 

247 

0 

.00000 

6.033E-02 

.10150 

.08799 

5.254E-12 

5.887E-12 

6.000E-02 

4.800E-01 

2.50E*01 

1 .40E*0l 

549 

211 

0 

.00000 

7.062E-02 

.09650 

.06456 

6.010E-12 

1 .679E-11 

7 -OOOE-02 

5.600E-01 

2 .5OE*01 

1 .40E*01 

497 

185 

0 

.00000 

8.055E-02 

.06850 

.07769 

1 .339E-U 

1 .282E-10 

8-000E-02 

6.400E-OI 

2 .50E*01 

1 .40E»01 

469 

164 

0 

.00000 

9.086E-02 

,08650 

.07630 

1 .935E'10 

7.669E-10 

9.000E-02 

7.200E-OI 

2.50E*01 

1 .40E*0l 

437 

148 

0 

.00000 

1 .007E-01 

.08600 

.07650 

1 .231E-09 

3.828E-09 

1 .OOOE-Ol 

8. OOOE-Ol 

INTEGRATION METHOD: 

SSFSFE 




















MAXIMUM 








NSTPR/ 

AVERAGE 



END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE CP 

- TIKE 

OVHD 

ERROR 

ERROR 

STEP SIZE 

STEP SIZE 

2.50E^01 

1 .40E*01 

1069 

495 

0 

.00000 

3.010E-02 

.16000 

.13632 

5.526E-12 

5.529E-12 

3.000E-02 

2.400E-0I 

2.50E*01 

1 .40E»01 

855 

371 

0 

.00000 

4.016E-02 

.14400 

.12541 

6.780E-I2 

6.789R.-12 

4.000E-02 

3.2Q0E-0I 

2.5OE^0l 

1 .40E^01 

721 

297 

0 

.00000 

5.017E-02 

.11200 

.09632 

7.145E-12 

7.221E-12 

5,O00E’02 

4. OOOE-Ol 

2.50E^0l 

1 40E^01 

621 

247 

0 

.00000 

6.0I1E-02 

.10400 

.09049 

5.254E-12 

5.8B7E-12 

6.000E-02 

4.600E-01 

2.50E»0l 

1 .40E»01 

549 

211 

0 

.00000 

7.0&2E-02 

.08850 

.07656 

6.010E-12 

1 .679E-U 

7 .OOOE-02 

5.600E-01 

2.50E^0l 

1 .40E«01 

497 

165 

0 

.00000 

8.055E-02 

.08000 

.06919 

1.339E-11 

1 .282E-10 

S.OOOE-02 

6.400E-01 

2.50E^0l 

1.40E^01 

469 

164 

0 

.00000 

9.086E-02 

.06550 

,07530 

1 .935E-10 

7.669E-I0 

9.000E-02 

7.200E-01 

2.50E^01 

1.40E«01 

437 

148 

0 

,00000 

1.007E-01 

.08050 

.07100 

1.231E-09 

3.828E-09 

1 .OOOE-Ol 

8. OOOE-Ol 


Table IV.I.AA 

COMPAR Summary of Statistics for the 
Euler Rigid Body Problem F*F(T,Y) 





8.S 


INTFCKATU'N MfcTHOl) S 7{S^ 


MAXIMUM 







nstpr; 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ABSERR 

RELERR 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SUE 

CP 

- TIME 

OVHD 

ERROR 

ERROR 

STEP SUE 

STEP SUE 

1 .OOE-04 

1 .OOE-12 

287 

19 

3 

.13636 

7.843E-01 


.03000 

.02456 

1 .912E-03 

1 .912E-03 

2 .OOOE-02 

1 .204E*00 

1 -OOE-06 

1 .OOE-12 

443 

31 

3 

.08824 

4.807E-01 


.05133 

.0429** 

4 .649E-06 

4-649E-06 

2 .OOOE-02 

6.405E-01 

I .OOE-08 

1 .OOE-12 

664 

50 

1 

.01961 

2.980E-0I 


.07700 

-06443 

6.091E-07 

6.091E-07 

2,Q00t-02 

4.531E-01 

1 .OOE-IO 

1 .OOE-12 

1106 

85 

0 

.00000 

1 .753E-01 


.10600 

.08505 

3.067E-09 

3 .067E-09 

2 .OOOE-02 

2 .791E-01 

I .OOE-U 

1 .OOE-13 

1457 

112 

0 

.00000 

1 -330E-01 


.15400 

.12641 

9.754E-11 

9.754E-11 

2 .OOOE-02 

2.114E-01 

1 .OOE-I2 

1 .OOE-13 

1925 

148 

0 

.00000 

1 .0U7E-01 


.19900 

.16254 

9.586E-12 

9.586E-12 

2.000E-02 

1 .584E-01 

1 .OOE-13 

1 .OOE-14 

2562 

197 

0 

.00000 

7.564E-02 


.26900 

.22048 

8.508E-13 

8.508E-13 

2 .OOOE-02 

1 .180E-OI 

integration METHOD: 

RK( 7)8 






















MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ABSERR 

RELERR 

npe 

NSTPA 

NSTPR 

NSTPT 

STEP SUE 

CP 

- TIKE 

OVHD 

ERROR 

ERROR 

STEP SUE 

STEP SUE 

1 .OOE-04 

1 .OOE-12 

282 

17 

5 

.22727 

8.765E-01 


.02240 

.01706 

1 .399E-03 

1 .399E-03 

2 .OOOE-02 

1 .192E*00 

1 .OOE-06 

1. OOE-12 

362 

25 

3 

.10714 

5.960E-01 


.02750 

.02064 

4.376E-05 

4.376E-05 

2.000E-02 

7 .770E-01 

1 .OOE-08 

1 .OOE-12 

631 

42 

7 

.14286 

3.548E-01 


.05200 

.04005 

1 .599E-07 

1.599E-07 

2.000E-02 

4.670E-01 

I .OOE’IO 

1 .OOE-12 

995 

70 

7 

.09091 

2.129E-OI 


.07950 

.06066 

1 .435E-09 

1 .435E-09 

2 .OOOE-02 

2.845E-01 

l.OOE-U 

1. OOE-13 

1270 

92 

6 

.06122 

1 .620E-01 


.09650 

.07245 

I.848E-10 

1 .848E-10 

2 .OOOE-02 

2.1971.-01 

1 .OOE-12 

1 .OOE-13 

1610 

120 

4 

.03226 

I .242E-0I 


.12500 

.09451 

2.551E-11 

2.551E-II 

2 .OOOE-02 

1 .710E-0I 

1 .OOE-13 

1 .OOE-14 

1952 

150 

0 

,00000 

9.934E-02 


.14500 

.10803 

3730E-12 

3.730E-12 

2 .OOOE-02 

1 .339E-01 

INTEGRATION METHOD: 

ODE 























MAXIMUM 








NSTPR/ 

AVERAGE 




END POINT 

GLOBAL 

MINIMUM 

MAXIMUM 

ABSERR 

RELERR 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SUE 

CP 

- TIME 

OVHD 

ERROR 

ERROR 

STEP SUE 

STEP SUE 

1 .OOE-04 

1 .OOE-12 

175 

87 

0 

.00000 

1 .713E-01 


.05650 

.05319 

7.061E-04 

6.772E-04 

2.039E-03 

2.610E-01 

l.OOE-06 

1 .Ov.E-12 

311 

155 

0 

.00000 

9.6UE-02 


.09500 

.08911 

1 .903E-05 

1 .897E-05 

2 .039E-04 

1 .044E-01 

1 OOE-08 

1. OOE-12 

437 

218 

0 

.00000 

6.835E-02 


. 1 3800 

.12972 

6.829E-08 

6.807E-08 

2.039E-05 

8.352E-02 

1 .OOE-10 

1 .OOE-12 

595 

297 

0 

.00000 

5.017E-02 


. 1 8800 

.17673 

1 .006E-09 

1 .OlOE-09 

2.039E-06 

6.68IE-02 

1 .OOE-n 

1 .OOE-13 

769 

384 

0 

.00000 

3.881E-02 


.23900 

.22444 

1 .654E-10 

9.693E-07 

6.448E-07 

4.226E-02 

1 .OOE-12 

1 .OOE-13 

865 

432 

0 

.00000 

3.449E-02 


.29400 

.27762 

4.833E-12 

9.196E-07 

2 .039E-07 

5.345E-02 

! .OOE-13 

1. OOE-14 

1133 

566 

0 

.00000 

2 .633E-02 


.37000 

.34854 

5.965E-12 

1 .454E-06 

6.448E-08 

3.380E-02 


INTCGRATION METHOD: EROCH2 


V 

r 


ABSERR RELERR 


NSTPR/ AVERAGE 

NFE HSTPA NSTPR HSTPT STEP SUE CP - TIME 


MAXIMUM 

END POINT GLOBAL MINIMUM MAXIMUM 
OVHD ERROR ERROR STEP SIZE STEP SIZE 


1 .OOC-04 

1 .OOE-12 

2 78 

U1 

0 .00000 

1 .057E-01 

.10800 

.10274 

1 .016E-03 

1 .OOE-06 

1 .OOE-12 

384 

199 

0 .00000 

7.488E-02 

.14000 

.13273 

3.674E-07 

l.OOE-08 

1 .OOE-12 

561 

284 

0 .00000 

5.247E-OZ 

.21700 

.20638 

2.767E-08 

1 .OOE-IO 

1 .OOE-12 

763 

390 

0 .00000 

3.821E-02 

.29500 

.28055 

7.537E-10 

1 .OOE-ll 

1.00E-. j 

798 

41C 

0 .00000 

3.582E-02 

.31400 

.29889 

6. 74 IE- It 

1 .OOE-12 

I .OOE-13 

817 

428 

0 .00000 

3.482F.-02 

.’MOO 

.29553 

1 .184E-U 

1. OOE-13 

1 .OOE-14 

1542 

783 

0 .uOOOO 

1 .903E-02 

.5b 300 

.55380 

7.319E-12 


1.017E-03 l,O00E-02 I .6O0E-OI 
S.B60E-07 6.7^0E-04 1 .600E-01 
2.767E-08 7.813E-OS 8.000E-02 
7.494E-10 9.766E>0b 8.000E-02 
6.726E-11 2.441E’06 4.000E-Q2 
9.176E-07 6.I04E-07 JOE-02 
9.I76E-07 3.0^2E-07 2.000E-02 


lltTECRATlON METHOD: KROCNl 


ABSERR 

RELERR 

NFE 

NSTPA 

NSTPR 

NSTPR/ 

MSrPT 

AVERAGE 
SlEP SUE CP 

- TIKE 

OVHD 

END POINT 

ERROR 

ii 

MINIMUM 
STEP SUE 

MAXIMUM 
STEP SIZE 

l.OOE-04 

l.OOE-12 

269 

138 

0 

.00000 

1.080E-01 

.12500 

.11991 

4 .864E-04 

4.849E-04 

1 .OOOE-02 

1 .600E-01 

1 .OOE-06 

1 .OCE-12 

391 

200 

0 

.00000 

7.451E-02 

'200 

.16460 

1 .830E-07 

1 .836E-07 

6.250E-04 

I .600E-01 

1 .OOE-08 

1 .oor-12 

567 

287 

0 

.00000 

5.192E-02 

..*5900 

.24826 

2 .296E-0B 

2 .296E-08 

7.813E-05 

8.000E-02 

1 ,OOE-i0 

l.OOE-12 

780 

399 

0 

.00000 

3.:35E-02 

.36000 

.3452? 

1 .254E-10 

1 .235E-10 

9.766E-06 

8.000E-02 

1. .E-ll 

1. OOE-13 

798 

416 

0 

.00000 

3.582E-02 

.37800 

.36289 

9.761E-12 

9.831E-12 

2.441C-06 

4.000E-02 

1 .OOE-12 

i. OOE-13 

IU26 

529 

0 

.00000 

2.817I-D2 ‘ 

.50200 

.48257 

1.305E-11 

9.176E-07 

6. 104E-07 

4.000E-02 

1. OOE-13 

1. OOE-14 

1542 

783 

0 

.00000 

1.903E-02 

.70200 

.67280 

4.7401-12 

9.176E-07 

3.052E-07 

2.000E-02 


Table IV. 1 . 4 b 


COMPAR Summary of Statistics for the 
Euler Rigid Body Problem F*F(T,Y,Y^^^) 





86 


i 


INThf.RATIOK MtTHOD ABFS 


ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPR/ 

NSTPT 

average 
STEP SIZE 

CP 

- TIME 

2 -30E*01 

1 .lOE^Ol 

1123 

496 

0 

.00000 

3.004C-02 


. 1 7000 

2.30E»01 

1 .IOE*01 

838 

372 

0 

.00000 

4 .006E-02 


.13800 

2.3OE*01 

1 .lOE^Ol 

740 

298 

0 

.00000 

3 .OOOE-02 


.10700 

2.30E*01 

1 .lOE+01 

631 

248 

0 

.OOJOO 

6.009E-02 


.10000 

2.30E*01 

1.10E*0l 

390 

212 

0 

.00000 

7.029E-02 


.09400 

2.30E^01 

1 .lOE+01 

3 38 

186 

0 

.00000 

8.011E-02 


.08350 

2 .30E*01 

1 .lOE+01 

307 

165 

0 

.00000 

9.031E-02 


.07730 

2.30E401 

1 -10E*01 

473 

149 

0 

.00000 

1 .OOOE-01 


.07600 

INTEGRATION METHOD: 

KSCFS 












NSTPR/ 

AVERAGE 



ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE 

CP 

- TIME 

2.30E*01 

1 .lOE^Ol 

1114 

496 

0 

.00000 

3.004E-02 


.13400 

2.30E*01 

1 .lOEoOl 

877 

372 

0 

.00000 

4.006E-02 


.11700 

2.30E*0l 

1 .lOE^Ol 

740 

298 

0 

.00000 

3 .OOOE-02 


.10700 

2.30E*01 

1 .tOE^OI 

640 

248 

0 

.00000 

6.009E-02 


.09200 

2.30E*0i 

1 -10E*01 

390 

712 

n 

.00000 

7 .029E-02 


.09500 

2.30C^)1 

1 .lOE^Ot 

338 

186 

0 

.00000 

8.01 lE-02 


.08800 

2.30E*0l 

1 .lOE^Ol 

496 

163 

0 

.00000 

9.0J1E-02 


.C80OO 

2.30E*01 

l.lOE^Ol 

464 

149 

0 

.00000 

1 .OOOE'Ol 


.08000 


INTEGRATlOfl KETHUD: SSFSBD 


NbTPR/ AVERAGE 


ITER 

ORDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE CP 

- TIME 

2.30E^01 

: -OOE^Ol 

1103 

496 

0 

.00000 

1 .004E-02 

.14403 

2.30E*0; 

1 .OOE^Ol 

863 

372 

0 

.00000 

4.006E-02 

. 1 2900 

2.30E«01 

1 .OOE^Ol 

727 

298 

0 

.00000 

S.OOOE-02 

.10350 

2.30E«01 

1 -OOE+Oi 

627 

248 

0 

.00000 

b.c: '1-02 

.08700 

2.30E*01 

1 .OOE^Ol 

363 

212 

0 

,00000 

7.029E-02 

.38800 

2.30E»01 

I .OOE^Ol 

323 

1P6 

0 

.00000 

8.01 lE-02 

.07930 

2.30E*0I 

1 .OOE^Ol 

491 

163 

0 

.00000 

9 .031E-02 

.07200 

2.30E»01 

1 .OOE^Ol 

449 

149 

0 

.OOOOU 

1 .OOOE-01 

.07200 


IMTEGRATIOH HETHOO i SSFSFE 


nstpr; aveiucc 


ITER 

CMDER 

NFE 

NSTPA 

NSTPR 

NSTPT 

STEP SIZE CP 

- T7ME 

2.30£»01 

1.00E«01 

1103 

496 

0 

.00000 

3.004E-02 

. 1 2800 

2.30E«01 

1 .OOE^Ol 

863 

372 

0 

.00000 

4.006E-02 

.10400 

2.30E^01 

1 .OOC^Ol 

727 

298 

0 

.00000 

5.000E-02 

.10230 

2.30E»01 

1 .O0E«01 

627 

248 

0 

.00000 

6.009E-02 

.08330 

2.30E«0I 

1 .OOE^Ol 

363 

212 

0 

.00000 

7.029E-02 

.07900 

2.30E*01 

1.00E«01 

523 

186 

0 

.03000 

S.OllE-02 

.07930 

2.30E«01 

t.OOE^Ol 

491 

163 

0 

.00000 

9.031E-02 

.07430 

2.30t»01 

1 .00E«01 

449 

149 

0 

.00000 

1 .OOOC-01 

.07000 


“AV IHI’M 

END POINT GLOBAL MINIMUM MAXIMUM 
OVHD ERROR ERROR STEP SIZE STEP SIZE 


.14869 9.061E-13 9.637E-13 3.000E-02 1 .8OOE-OI 
.12113 3.636E-12 6.273E-12 4.000E-02 2.400E-01 
.09299 6.861E-11 7.888E-11 5 .OOOE-02 3.000E-01 
.08767 3.289E-10 6.013E*10 6.000E-02 3.600E-01 
.08283 3.045E-09 3.200E-09 7.000E-02 4.200E-01 
.07331 1.O78E-O0 1.474E-08 0.OOOE-O2 4.800E-01 
.06790 8.679E-07 7.415E-07 9.000E-02 5.400E-01 
.06700 1.359E-04 8.401E-05 l.OOOE-OI 6.000E-01 


MAXIMUM 

END POINT GLOBAL MINIMUM MAXIMUM 
OVHD ERROR ERROR STEP SIZE STEP SIZE 


.13290 6.672E-12 6 710E-12 3 OOOE-02 I .800E-01 
.10039 7.234E-I2 7.061E-12 4.000E-02 2.400E-01 
.0‘J299 1.089E-10 1.069E-10 3.000E-02 3.000E-01 
.0 ‘»8 1.304E-09 1.300E-09 6.000E-02 3.600E-01 
.08. j 1.019E-08 1.029E‘08 7.000E-02 4.200F.-01 
.07781 5.972E-08 6.000E-08 8.000E-02 4.800E-01 
.07861 2.751E-07 2.78U-07 9.000E-02 5.400E-01 
.07121 1.040E-06 I.115E'06 1 .OOOE-01 6.000E-01 


MAXIMUM 

END POINT GLOBAL MINIMUM MAXIMUM 
OVHD ERROR ERROR STEP SIZE STEP SIZE 


.12311 2.223E-12 2.285E-12 3.000r-02 1 .830E-0I 
.11262 1.854E-t2 2.138E-12 4.000E-02 2.400E-01 
.06973 1.869E<12 1 .346E’ll S .QOOE-02 3.000E-01 
.07313 1.075E-10 I-272E-10 6.000E-02 3.600E-01 
.07730 1.2S0E‘09 1.318E-09 7.000E-02 4.2OOE-0I 
.06960 8.47IE-09 8.770E-09 8.000E-02 4.800E-01 
.06270 4.220E-08 4.368E-08 9.000F-02 3.400E-01 
.06330 1.710E-07 2.796E-07 1 .OOOE-01 6 .OOOR-Ol 


MAXIMim 

END POINT GLOBAL MINIMUM MAXIMUM 
OVHD ERROR ERROR STEP SIZE STEP SIZE 


.10711 2.223E-12 2.283E'12 3.000E-02 1.800E-01 
.08767 1.634E-12 2.138E-12 4.000E-02 2.400E-01 
.08873 1.869E'12 1.346E*11 3.000E>02 3.000E-01 
.07163 1.073E-10 1.272E'lO 6.000E-02 3.600E-01 
.06830 1.230E-09 1.318E-09 7.000E-02 4.200E-01 
.06960 8.471E-09 6.770E-09 S.000E<02 4.600E'01 
.06320 4.220E>08 4.368E'08 9.000E-02 3.4COE-01 
.06130 1.710E-07 2.796E'07 1 .OOOE-Ol 6.000E-01 


Table IV.1*4B 

COMPAR Summary of Statistics for the 
Euler Rigid Body Problem F-F(T,Y,Y^^^) 
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IV. 1.5 Rest ricted Thrco-Pody Problem 


The restricted three-body problem in two dimensions may be 
described by the set of first-order differential equations 




^ = Y 4- Y — — ! 

dt *1 3 3 

•"l '‘2 


dx^ x^Cl-VJ) Ux_ 

"dr ' "^*3 * ^2 2 T 

'"l '‘2 


2 2 2 2 2 5 

where = (x^+p) +X2 t - (x^>p-1) ♦y ' and U is the mass ratio. 

By selecting the initial conditions x^(0)=i.2, X2(0)=0.0, x,(0)=0.0 

and Xjj(O) = -1.04935750983031990726 and the mass ratio as 

y = 1.0/82,45, the period of the orbit becomes 

T = 6. I92I6933131963970674. Although the state at T is equal to the 

initial conditions, i.e., x(T) = x(0), there is no analytic solution in 

the interval 0 < t < T. The reference solution for this problem was 

generated by using ODE with an absolute error tolerance of 10“^^ and a 

-12 

relative error tolerance of 10 . The reference orbit is shown in 

Figure IV. 2. 5a. 


Several computer runs were made in an attempt to determine 
the optimum order for the fixed-mesh/fixed-order integrators. However, 
to obtain maximum global errors of less than 1.0, the range of possible 
stepsizes had to be restricted to less than 0.001. For a stepsize of 
0.001, the multistep integrators require approximately 12,500 function 
evaluations, while the most function evaluations required by any of the 


Jf 








- 


Efficiency Curves for the Restricted Problem of 

Three Bodies 
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INTKllKATlON METHOD' 

ABSEKK RELERt 

SFE 

7(8) 

NSTFA 

NSTPR 

NSTFR/ 

NSTPT 

AVERAGE 
STEP SUE 

CP - TIME 

OVHD 

END POINT 

ERROR 

MAXI.iUM 

CLORAL 

ERROR 

MINIMUM 
STEP SUE 

MAX I 

STEP SUE 

1 .OOE-02 

1 OOE-12 

104 

6 

2 

.25000 

1 .950E-01 

0.00000 

0.00000 

METHOD PARED TO REACH PINAL TIME 

1 .OOE-04 

l.OOE-12 

833 

44 

20 

.31250 

1 .407E-01 

.09500 
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6.725C*00 

1 .OOOE-03 

4.9671-01 

1 .OOE-06 

1 .OOE-12 
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55 

21 

.27632 
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76 

22 
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.14300 

.11927 
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1.786E-0) 

1 .OOOE-03 
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1612 
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12 
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.16126 
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\ .OOE-ll 
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5 
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.18673 
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1 .0741-01 
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MINIMUM 

MAXIMUM 

ARSERR 

RELERR 

N»E 

HSTPA 
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NSTPT 

STEP SUE 

CP - TIMR 
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ERROR 

ERROR 

STEP SUE 
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97 

6 

1 
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47 

13 
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1 185 

68 

25 
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5. 8531-01 
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42 
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3.5321-01 
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63 

.26250 
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.23600 

.17909 
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7 .610E-04 
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1 .OOE-ll 

1 .OOE-13 

3917 

2)1 

76 

.24756 
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.29700 

.22411 
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2.395E-07 
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1.4781-01 
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MINIMUM 

MAXIMUM 
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NSTPT 

STEP SUE 
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OVHD 

ERROH 

ERROR 

STEP SUE 

STIP SUl 
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1 .OOE-12 

346 
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0 
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3.846E-02 
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.08006 

1 .76OE^O0 

6.594E«00 

1 .073C-U3 

2.7481-01 

1 .OOE-04 

1 OOE-12 

516 
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0 
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2.517E-02 

.15200 

.14240 

6.421E-03 
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5.634E-04 

1 .<»421-0l 

1 .OOE-04 

1 .OOE-12 

814 
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0 
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1 .556E-02 

.24200 

.22685 

5.417E-05 

3.726E-02 

1 .718E-04 

1 .6821-01 

1 .OOE-08 

1 .OOE-12 

1142 

561 

0 
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1 .104E-02 

.37900 

.35775 

6.435E-07 

5.497E-04 

1.718C-05 

7.4771-02 

1 .OOE-10 

1 .OOE-12 

i671 

826 

0 

.00000 

7.497E-C3 

.54200 

.51090 

8.816E-09 

7.607C-06 

1.7181-06 

7.0641-02 

1 .OOE-ll 

1 OOE-1) 

1939 

962 

0 

.00000 

6.437E-0) 

.61200 

.57592 

3.040E-10 

1.1511-06 

5.43IE-07 

5.2171-02 

IHTECRATIOR METHOD: 
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NSTPR/ 

AVERAGE 



CHD POINT 

CLORAL 

MINIMUM 

MAXIMUM 

ARbERI 

■.^LERR 

HFE 

HSTFA 

NSTPR 

NSTPT 

STEP 

CP - TIME 

OVHD 

ERROH 

ERROR 

STPP SIZE 

STIP SUt 

i .OOE-02 

1 .OOE-12 

456 

223 

0 
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2.777E-02 

.14300 

.13451 

3.882E-02 

4.3821*00 

1 .OCOE-03 

2. 5601-01 

1 .OOE-04 

1 .OOE-12 
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348 

0 
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1.7791-02 

.22200 

.20897 

2.170E 04 

2.498E-01 

5.000E-04 

1.2801 1 

1 .OOE-04 

1 OOE-12 

i 105 
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0 
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1.122E-02 

.35900 

.3)844 

5.165E-06 
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2.5001-04 
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1 OOE-12 
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0 
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.16600 
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1 .0001-03 
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l.OOE-12 
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.24900 

•01 
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I.OOE-04 
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109) 
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l.OOE-12 
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.00000 
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4.235E-08 
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1 OOE-10 

1 OOE-12 
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0 
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.70771 

I.700E-10 

T.275E-07 
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1. OOE-ll 
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.90700 

.86176 
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).:oot-o2 


Taulcf IV. 1.5 

COMPAR Summary cf Statistics for the 
Restricted Problem •'[ Three Bodies 


/ 
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var iable-st«?p integrators was 30"^? and, In general, less than 2000. 
Since the pcrfornances of the fixed-nesh/fixed-order inte.-^rators are 
inferior to the variable-step integrators for this problem, the effi- 
ciency curves of the fixed-mesh/fixed-order integrators are r.Jt 
presented. Figures IV. 1.5b through IV.1.5e show the efficiency curve 
for ODE, PK(7)8, RKN7(8), KROGHI and KR0GH2 for an integration of one 
period. The results of the variable-step integrators are summarized in 
Table IV. 1.5. 

IV. 1.6 Conunentf on the COMPAR Results 

In general., the efficiency curves indicating the maximum glo- 
bal errors for ODE, KROGH i and HR0GH2 increase dramatically for the 
more strict absolute and relative error tolerances. The maximum global 
error for high tolerances occurs within the first few .nteps in the 
integration by ODE, KROGHI or KR0GH2, i.e., during Ihe bootstrapping 
starting procedure of these integrators when the order of the integra- 
tion formulas is being increased. 

The maximum stipsizes listed in Tables IV. 1.1 th.'ough TV.1.M 
for the fixed-mesh raultlstep integrators result from the procedure COM- 
PAR used t'' calculate the maximum steps! ze and the starting algorithm 
used by the fixed-mesh integrator. As a result, the maximun stepsize 
listed for the fixed-mesh integrators is toe difference in the value of 
the ind^'pensent variable at the initial conditions and at the last node 
point in the starting procedure. 
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IV. 2 Comparison of Integration Errors for a Typical Satellite Problem 

*^e major focus of this report is the evaluation of the 
fixed-mesh/ fixed-order multistep integrators for solving the satellite 
problem. The evaluation is carried out by comparing the results of the 
fixed-step integrators with the results of the variable-step integra- 
• -rs for two different force models: a two-body model and an eleventh 
degree and order geopotential model. The comparisons are carried out 
for an integration interval of 30 days. In order to reduce the number 
of parameters that could be varied during the evaluation, only one set 
of initial conditions was examined. This set of initial conditions 
corresponds to the Lageos satellite. The initial conditions in rec- 
tangular coordinates are 

• 

X = -9123000.0 m x = 2708.0 m/sec 

y = -1894000.0 m y = 3145.0 m/sec 

z s -7822000.0 m z ~ -3930.0 m/sec 

and the osculating set of Keplerian orbital elements are 

a s 1.226538 x 10^ m 12 * .5092370 radians 

e = .003530571 0) s -1.6593 11 racians 

i s 1.918821 radians H s 5.553645 radians 





The comparisons for each force model are carried out in rec- 
tangular coordinates, orbital elements and radial-transverse-normal 


r 


V 
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(RTN) coordinates relative to the reference orbits. The errors are 
computed by subtracting the integrated solution from the reference 
solution. 


The RTN coordinates are defined relative to the position and 
velocity vectors of the reference orbit. The radial unit vector Is 
defined to be in the direction of the positon vector. The normal unit 
vector is both parallel to and in the direction of the angular momentum 
vector. The transverse unit vector is defined to be orthogonal to the 
normal and radial unit vectors in the approximate direction of the 
velocity vector. 

An examplr of a set of comparison plots are shown in Figures 

IV. 2.0a through IV. 2. Or. The differences illustrated in Figures 

IV. 2.0a through 2. Or are the differences of integrating a spherical 

eleventh degree and order model with KSGFS using a stepsize of 300 

seconds and an order of 14 and integrating the same model with a 

double-precision RK(7)8 using absolute and relative error tolerances of 
—1 8 

10 . To compare the six components of each of the three coordinate 

systems would be cumbersome and probably ambiguous. Since it has been 
noted that a major portion of the integration error is re.'lected in the 
transverse error, the transverse error results were chosen as the basis 


for the comparisons. 
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IV. 2.1 Two - Body Force Model 

The reference orbit for the two-body model was generated by 

using the analytical solution of the two-body problem using orbital 

elements. The Keplerian elements a, e, i» ^ and u are constant for 

the two-body, point-mass problem. The mean anomaly, M , was calculated 

from the equation M s M nt where n is the mean motion and M is 

o o 

the value of the mean anomaly at t s 0. 


The first-order ordinary differential equations for this 
problem were 



and the second-order differential equations are 



. 2 2 2 2 
where r * + Xg ♦ . 



• i* f 
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A plot of the transverse errors of the variable-step integra- 
tors ODE, RK(7)0, RKN7(8), KR0GH1 and KR0GH2 are shown in Figure 

IV. 2. la. An absolute error tolerance of and a relative error 

-12 

tolerance of 10 was used for ODE, RK(7)8 and RKN7(8), and for 

_g 

KR0GH2, an absolute error tolerance of 10 and relative error toler- 
ance of 10“'^ was used. It should be noted that, using the most 
stringent tolerances allowed, the transverse error for KR0GH1 is 
approximately two orders of magnitude larger than the results from any 
of the other variable-step integrators. using the most stringent 
tolerances allowed by KR0GH1. 

To determine the behavior of the transverse error from the 
fixed-step integrators, they were applied with a variety of stepsizes 
and orders to integrate this problem. Table IV. 2.1 is a summary of the 
results and lists the maximum transverse error encountered for each 
integrator for each combination of order and stepsize. Figure IV. 2. 1b 
shows the transverse error curves of ABFS for an integration order of 9 
and stepsizes of 50.0, 75.0, 100.0 and 150.0 seconds, while Figure 
IV. 2. 1c shows the transverse error curves of KSGFS for an integration 
order of 13 and stepsizes of 100.0, 200.0 and 250.0 seconds. 

IV. 2. 2 Spherical Eleventh Degree and Order Force Model 

This force model does not have an analytical solution and is 


.V 




characterized by the first-order ordinary differential equations 



\sl 




Figure IV. 2.1 a 

Transverse position error of the variable 
step and variable mesh int«igrators for 
the two body force model with an analytic 
t%ro body solution. 




Table IV. 2.1 


lOA 


Maximum Transverse Errors (m) 
Two Body Force Model 
'.hirty Day Integration Interval 


AaFS 

STEPSIZE/ORDER 

7 

8 

9 


10 

50 sec 

- .221 

-.097 

-.034 

- 

.048 

100 sec 

-16.90 

4.20 

-.005 

- 

.048 

150 sec 

- 642. 

159. 

3.65 


.767 

200 sec 

-8449. 

2035. 

85.91 


16.08 

KSG 

STEPSIZE/ORDER 

11 

12 

13 


14 

100 sec 

- .038 

-.041 

-.037 

- 

.059 

200 sec 

- .004 

-.057 

-.055 

- 

.043 

300 sec 

7.70 

0 
• 

1 

-.224 


.052 

400 sec 

260. 

11.59 

-7.64 

- 

1.55 

SSFSFE 

STEPSIZE/ORDER 

11 

12 

13 


14 

100 sec 

- .053 

-.035 

-.041 

- 

.024 

200 sec 

- .062 

-.045 

-.050 

- 

.050 

300 sec 

- .319 

-.077 

-.030 

- 

.051 

400 sec 

- 8.83 

-2.80 

.094 


.104 

SSFSBD 

STEPSIZE/ORDER 

11 

12 

13 


14 

100 sec 

- .064 

-.046 

-.040 

- 

.047 

200 sec 

- .042 

-.046 

-.043 

- 

.042 

300 sec 

- .088 

0 

e 

1 

-.042 

•• 

.046 

400 sec 

- 3.01 

-.206 

•f.108 


e 


0 
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dx^ 

“dt ' ’'m ■ " 


dXg 

dT = ’'5 = y 


ox. 


r 

px 


T * 8l 


*5 ^ 


♦ 8^ 


dx, 

“dt" ^ *6 = ^ 


Vix, 


*6 “ 


T* 83 


or the second-order ordinary differential equations 


d^, 


dt 




dt 


d^x- 


dt 


WXi 


r ♦ 8, 


Ux 

r 


2 

3* «2 


Ox. 


J" *3 


where r = x.j + Xg + x^ and g.^, and g^ are components of which 

is the acceleration vector, g , due to a nonspherical earth. The 

reference orbit for this model was generated by using a double- 

—1 8 

precision version of RK(7)8 with an absolute error tolerance of 10 

—1 8 

and a relative error tolerance of 10 . There were 927,317 function 

evaluations made during the generation of the reference orbit for an 
orbital arc of 2,700,000 seconds. To verify that the tolerances used 
in generating the reference orbit were not too strict, the reference 
orbit was compared to another integration by the double-precision 
RK(7)8 that used an absolute error tolerance of 10”^^ and a relative 


tolerance of 10 


-15 


This second integration required 376,897 function 
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evaluations to integrate an orbital arc of 2,500,000 seconds. The com- 
parison of these two trajectories was done in orbital elements and RTN 
components which are shown in Figures IV. 2. 2a through IV. 2. 21. The 
differencing of these two comparison trajectories was done in single- 
precision, and Figures 4.2.2a through 4.2.2f indicate differences in 
the two orbits at the level of the roundoff error of the computer, 
except for a small secular difference in the mean anomaly. The differ- 
ences in the two orbits are magnified when the radial, transverse and 
normal differences are examined. Despite the growth of the transverse 
position difference between the orbits, the difference is still rela- 
tively small. Thus, from an accuracy and roundoff error point of view, 

—1 8 

restricting the tolerances to 10~ is not inappropriate. 

The fixed-step integrators ABFS, KSGFS, SSFSBD and SSFSFE use 
the PECE algorithm as described in Chapters II and III. It has been 
noted that if the value of the state does not change significantly from 
the predicted value to the corrected value, then the value of the per- 
turbing function, g , may vary only slightly from the prediction step 
to the correction step. The calculation of g generally requires many 
more arithmetic operations than the two-body terms. The calculation of 
g for each calculation of the two-body terms is referred to as a "full 
evaluation." However, if g does not change significantly between the 
prediction step and correction step, then g need only be calculated 
during the prediction step. The calculation of g only once per 
integration step is referred to as a "partial evaluation." 



r 
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Figures IV. 2. 2 a,b,c 

Differences in Semimajor 
Axis, Eccentricity, and 
Inclination for Spherical 
Eleventh Degree and Order 
Model 

Double Precision RK(7)8 using 10~ 
tolerances vs Double Precision 
RK(7)8 using lO"^*^ tolerances 
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Figures IV. 2. 2 d,e,£ 

Differences In Ascending 
Node, Arguement of Perigee, 
and Mean Anomaly for Spher- 
ical Eleventh Degree and 
Order Model 

Double Precision RK(7)8 using lO"^^ 
tolerances vs Double Precision 
RK(7)8 using 10 tolerances 
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Figures IV. 2. 2 g.h.i 

Radial, Transverse, and 
Normal Position Differences 
for a Spherical Eleventh 
Degree and Order Model 

Double Precision RK(7)8 using 10 ^ 
tolerances vs Double Precision 
RK(7)8 using lO" tolerances 
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Figures IV. 2. 2 j,k,l 

Radial, Transverse, and 
Normal Velocity Differences 
for a Spherical Eleventh 
Degree and Order Model 

Double Precision RK(7)8 using 10 ^ 
tolerances vs Double Precision 
RK(7)8 using 10 tolerances 
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Figure IV. 2. 2m is a plot of the 

transverse errors 

of the 

variable-step 

integrators, ODE, RK(7)8, 

RKN7(8), 

KR0GH1 and 

KR0GH2. 

The tolerances 

used for each integrator were 



INTEGRATOR 

ABSOLUTE ERROR TOLERANCE 

RELATIVE 

ERROR TOLERANCE 

ODE 

10-1° 

10-12 



KR0GH1 

10"^ 

10-^ 



KR0GH2 

10"^ 

10-11 



RK(7)8 

10-1° 

10-12 



RKN7(8) 

10-1° 

10-12 . 




As in Section IV.2.1, the fixed-step integrators were used 
with a variety of stepsizes and orders to integrate this problem. 
Tables IV.2.2A and IV.2.2B are the summaries of the results of each 
integrator for each pair of order and stepsize. Table IV.2.2A are the 
results when full evaluations are made, and IV.2.2B are the results 
when partial evaluations are made. 

To illustrate the behavior of the transverse error for each 
of the fixed-step integrators, Figures IV.2.2n through IV.2.2q are the 
plots of the transverse errors for ABFS, KSGFS, SSFSBD and SSFSFE with 
integration orders of 9, 13, 13 and 12, respectively, with a variety of 
stepsizes and with full evaluations of g . 






integrators solving the spherical eleventh order and 
degree problem: 1) with KROGHl and 2) without KROGHl, 





Table IV.2.?A 
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Maxinum Transvcrrc Trror (m) 

Spherical Eleventh Degree and Order Ccopotentlal 
Thirty Day Integration Interval 
Full Evaluation of Accelerations 


ABFS 

STEPS I ZE/ ORDER 

7 

8 

9 

10 

50 sec 

- .032 

-.092 

-.026 

-.039 

100 sec 

- 16.9 

4.17 

-.019 

-.053 

150 sec 

- 639. 

158. 

3.60 

-.750 

200 sec 

-8421. 

2026. 

85.0 

-15.7 

KSG 

STEPSIZE/ORDER 

11 

12 

13 

14 

100 sec 

- .025 

-.056 

-.044 

-.046 

200 sec 

- .019 

-.019 

-.016 

-.031 

300 sec 

5.66 

.170 

2.60 

-.683 

400 sec 

231. 

-9.85 

14.9 

31.3 

SSFSFE 

STEPSIZE/ORDER 

11 

12 

13 

14 

100 sec 

- .052 

-.053 

-.046 

-.040 

200 sec 

- .045 

-.039 

-.035 

-.054 

300 sec 

- .191 

.104 

.032 

-.097 

400 sec 

- 6.51 

1.72 

3.34 

-2.77 

1 

SSFSBD 

: :£PSIZE/0RDER 

11 

12 

13 

14 

100 sec 

- .044 

-.040 

-.055 

-.039 

200 sec 

- ,042 

-.040 

-.045 

-.036 

300 sec 

.053 

.068 

.009 

.024 

400 sec 

- .480 

1.68 

1.13 

* 


¥ 






Table IV.2.2B 
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Maximum Transverse Error (m) 

Gjilierical i^levcnth Degree and Order Geopotential 
Thirty Day Integration Interval 
Partial Evaluation of Accelerations 


ABFS 

STEPSIZE/ORBER 

7 

8 

9 

10 

SO sec 

- .025 

-.090 

-.035 

-.045 

100 sec 

- 16.9 

4.17 

-.004 

-.066 

150 sec 

- 641. 

158. 

3.65 

-.748 

200 sec 

-8431. 

2027. 

85.9 

-15.7 

KSG 

STEPSIZE/ORDER 

11 

12 

13 

14 

100 sec 

- .040 

-.058 

-.025 

-.031 

200 sec 

- .029 

-.038 

-.019 

-.043 

300 sec 

5.05 

.360 

2.74 

-.707 

400 sec 

196. 

-13.7 

23.7 

35.0 

SSFSFE 

STEPSIZE/ORDER 

11 

12 

13 

14 

100 sec 

- .060 

-.046 

-.040 

-.045 

200 sec 

- .042 

-.043 

-.041 

-.047 

300 sec 

- .959 

.361 

.207 

-.140 

400 sec 

- 49.9 

-4.34 

14.38 

2.05 

SSFSBD 

STEPSIZE/ORDER 

11 

12 

13 

14 

100 sec 

- .054 

-.056 

-.050 

-.039 

200 sec 

1 

o 

0 

• 

1 

-.045 

00 

0 

e 

1 

300 sec 

- .706 

.320 

.165 

.014 

400 sec 

- 42.7 

-4.01 

11.9 

* 
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Figure IV. 2, 2 n 

Comparison of transverse errors from solving 
the spherical eleventh degree and order 
model with ABFS (order ■ 9) using various 
stepsizes. 


4 »' 


* 



U: 



Figure IV. 2, 2 o 

Omparison of transverse errors from solving 
the spherical eleventh degree and order 
model with KSCFS (order * 13) using various 
stepslzes. 
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Figure IV, 2, 2 p 

Comparison of transverse errors from solving 
the spherical eleventh degree and order 
model with SSFSBD (order • 13) using 
various stepsizes. 
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Figure IV. 2. 2 q 

Comparison of transverse errors from solving 
the spherical eleventh degree and order 
model with SSFSFE (order • 12) using 
various stepsizes. 




4 



CHAPTER V 


Conclusions 

Results from COMPAR give an indication of the relative advan- 
tages and disadvantages of the three groups of integrators: 1) 
variable-step Runge-Kutta integrators RK(7)8 and RKN7(8); 2) variable- 
mesh/variable-order multistep integrators ODE, KR0GH1 and KR0GH2; and 
3) fixed-mesh/fixed-order multistep integrators ABFS, KSGFS, SSFSBD and 

SSFSFE. 

The harmonic oscillator problem of Section IV. 1.1 and the 
circular two-body problem of Section IV. 1.2 have the same analytic 
solution but have linear and nonlinear differential equations, respec- 
tively. However, the relative performance of the three groups of 
integrators is the same in each problem. The fixed-mesh multistep 
integrators are the most efficient with respect to the number of func- 
tion evaluations required and central processor time used. The 
variable-mesh multistep integrators are competitive with the Runge- 
Kutta integrators in central processor time used and are more efficient 
than the number of function evaluations required. 

The elliptic two-body problem of Section IV. 1.3 illustrates 
the major differences between the variable-mesh multistep integrators 
and the single-step Runge Kutta integrators. The variable-mesh mul- 
tistep Integrators require fewer function evaluations than the 
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variable-step Runge-Kutta integrators to achieve a certain accuracy; 
however, the variable-nesh multistep integrators require more central 
processor time. Thus, the variable-mesh multistep integrators require 
more overhead per function evaluation than the variable-step Runge- 
Kutta integrators. If the function evaluations of a differential equa- 
tion are relatively inexpensive in central processor time, the Runge- 
Kutta integrators would have a distinct advantage; whereas, if the 
function evaluations are relatively expensive, the variable-mesh mul- 
tistep integrators would have an advantage. The performance of the 
fixed-mesh multistep integrators lie between the performances of the 
other two groups of integrators. The fixed-mesh multistep integrators 
reduce the overhead associated with the variable-mesh integrators while 
still requiring fewer function evaluations than the variable-step 
Runge-Kutta integrators to achieve a certain accuracy. 

The relative results of the three groups of integrators for 
the Euler rigid-body problem of Section IV. 1.4 are similar to the 
results for the eccentric two-body problem. The variable-step Runge- 
Kutta integrators require the least amount of overhead per function 
evaluation, while the variable-mesh multistep integrators require the 
fewest number of function evaluations to achieve a certain accuracy. 
Again, the fixed-mesh multistep integrators require less overhead per 
step than the variable-mesh multistep integrators and fewer function 
evaluations than the variable-step Runge-Kutta integrators. A point of 
interest with this problem is the relative performance of the Class 
I/fixed-mesh multistep integrator ABFS with that of the Class 



Il/fixed-mesh multistep methods. When the derivatives are a function 
of the states x^, x^ and x^ only, ABFS performs mo. e efficiently. How- 
ever, when the derivatives are a function of all six elements, the 
Class II methods perform as efficiently as ABFS. 

The performance of numerical integration algorithms to solve 
the satellite problem depends on many parameters, e.g., the distance of 
the satellite from the primary, the eccentricity of the satellite 
orbit, the intricacy of the model used to represent the primary. How- 
ever, by examining one scenario of the satellite problem, the charac- 
teristics of the fixed-mesh multistep integrators may be illustrated. 

Figures IV. 2. la and IV. 2. 2m are comparisons of the transverse 
errors for the variable-step and variable-mesh integrators. Table V.1 
is a summary of the approximate number of function evaluations required 
by each integrator. Figures IV. 2. la and IV. 2. 2m show that ODE and 
KR0GH2 are more accurate than RK(7)8 and RKN7(8), while Table V.l shows 
that ODE and KR0GH2 also required significantly fewer function evalua- 
tions. 


Tables IV. 2.1, IV.2.2A and IV.2.2B give the maximum 
transverse errers of the fixed-mesh multistep integrators for a variety 
of stepsizes and orders. These tables point out the characteristics of 
the formulations. First, it is noted that for stepsize less than 300 
seconds, the use of "partial evaluations" of the derivatives does not 
significan'cly alter the transverse errors. The second point is that 
the Class II formulations allow larger stepsizes than the Class I 
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Table V.l 

Approx invite Number of Function bv.iluatlons 
Required by F,ach Integrator for an Integration Interval 

of Thirty Days* 


INTEGRATOR; 

TWO BODY MODEL 

ELEVENTH DEGREE 
AND ORDER MODEL 

RK(7)8 

168,000 

168,000 

RKN7(8) 

770,000 

271,000 

ODE 

44,800 

53,900 

KR0GH2 

32,500 

58.900 

KROGHl 

30,300 

32,800 


FIXED MESH MULTISTEP INTEGRATORS WITH STEPS IZES OF: 


75 sec 

69,000 

69,000 

100 sec 

52,000 

52,000 

125 sec 

41,600 

41,600 

200 sec 

26,000 

26,000 

250 sec 

20,900 

20,900 

300 sec 

17,400 

17,400 


The tolerances used by the variable step integrators are given 
in Sections IV. 2.1 and IV. 2. 2. 
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formulation to achieve comparable accuracies. Finally, except for one 
combination of stepsize and order, there is no appreciable difference 
between the second-sum formulation using back differences and the 
second-sum formulation using function evaluations. 

The results of Section IV. 2 illustrate the advantages of 
fixed-mesh multistep integrators for solving the satellite problem. 
The fixed-mesh multistep integrators are capable of being as accurate 
as the variable-mesh multistep integrators while requiring fewer func- 
tion evaluations than the variable-mesh integrators. As noted above, 
the variable-mesh multistep integrators are more efficient than the 
variable-step Runge-Kutta integrators for solving the satellite prob- 
lem. 


The results of Section IV. 1 illustrate that the fixed-mesh 
multistep methods require less computer time overhead than the 
variable-mesh integrators. By coupling the results of Sections IV. 1 
and IV. 2, the fixed-mesh multistep integrators are shown to be an 
attractive tool for use in the satellite problem. 



APPENDIX A 


COEFFICIENTS FOR THE GENEKAL FORMULATION 


Equations (2.7) and (2.8) are the general formulation basic 

equations and are used to approximate the state at t given the state 

n+r 

and back differences at t^ . The coefficients, a and B • to be 
used in (2.7) and (2.8) are determined by the following algorithm:* 


1. Define h, = t - t 
I n+r n 



where h is the integration stepsize 


2. Calculate the matrix of coefficients • where 


for k = 1 


‘ q 


for k > 1 


®k,q ” ^k®k-1,q *^k®k-1,q+1 


where k = 1 , . . . ,i+l 
q - 1 , . . . ,k 


i = number of coefficients required in (2.7) and (2.8) 


= h + n - 2 

k s 


*This flxed-mesh coefficient algorithm Is determined from I'e variable- 
mesh coefficient algorithm given by Shamplne and Gordon (1975). 
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y = - U) 
k k ^k 


=rn 


1 

k 's 


3. The coefficients in (II. 7) and (II. 8) are given by 




“‘j.r = 82 ,j 





APPENDIX B 


COEFFICIENTS FOR THE SECOND-SUM FORMULATION 

Equations (2.24) and (2.25) are the second-sum formulation 

basic equations and are used to approximate the state at given 

the first and second sums and back differences at t . The coeffi- 

n 

cients, a and b , to be used in (2.24) and (2.25) are determined 
from the following recursive relationships:* 


b 


j-l.s 


J , 

1 y y 

k=o “ 


j=i 


i 


a 


J-2.S 


J n 

1 y y 

kfo k 'j-k.s 


j=2 i+1 


where i is the number of coefficients required in (2.24) and (2.25) 
and 


f 



k-1 

2 — • — y 

qfo •‘-‘I"’ ‘I 






1 

tn 


'in-1 ,s 


for m > 0 


*Thls is the algorithm for the second sum of coefficients as presented 
by Spier (1971). 
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where h is the integration stepsize 



APPENDIX C 


MODIFIED, GENERAL FORMULATION, PECE ALGORITHM 


Assuming the state (y^, and the i back differences of 
the function f(t,y,y^^^) are known at t^ , then the modified, general 
formulation, PECE algorithm that is used to advance the solution from 
tn to t^^^ is given by: 

1. Predict at t^^^ using an i^*^ order formula 


"n+1 


y -f hy^^^ + h^ 2 a f 

n ^n j,1 n 


Vi 




2. Evaluate the 
'»n.f "il!’ 


function f with the predicted solution 


n+1 


= f'Vr "ili’ 


3. Form the modified back differences 


s 2 f . ks1 i 

J=k " 


7^ fP 1 = fP , - d, 
n+1 n+l 1 


r 
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ni 


*1. Correct the solution using an (i-*-1) order formula 


y . = p , h^a. 1 1 1 

n+1 n+1 i+1,1 n+1 


+ h 6. . - fP , 

*^n+1 1+1,1 n+1 

5. Evaluate the function with the corrected solution 


^n+1 = ^^^n+r ^n+r ^n+1^ 


6. Advance the back differences from t to t , 

n n+1 


f i - d, 

0 n+1 1 


‘'n.l - <’<, * “o = ■’=' ‘ 



APPENDIX D 


MODIFIED, SECOND-SUM FORMULATION . PECE ALGORITHM 


-2 -1 

Assuming the second and first sums ( f , f ) and the 

n n 

1 1 ) 

i back differences of the function f(t,y,y' ) are known at t , 

n 

then the modified, second-sum formulation, PFCE' algorithm that is used 

to advance the solution from t tc t . is given by: 

n n>1 ^ ^ 

Predict (y, at t . using an i^^ order formula 


0„.1 - 


vJ"'r ) 




Evaluate the function 

<Pn.r 


with the predicted solution 


^n+1 ■ Pn>1* 


3. Form the modified back differences 


2 ; k-l i 


'‘Cl = Cl - "i 
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4. 'i->rrpnt the solution using an (i + 1) order formula 


.2 , i p 

y = p , ♦ h a. > , 

nf1 n+1 i,1 n*1 


K yieP 

Vi = Vi * ^ ^.1 ^ Vi 


5. Evaluate the function with the corrected solution 


^n+1 ' ^^^n+r ^n-t-T ^ni-1^ 


6, Advance the back differences from t to t - 

n n+1 


''o = Vi - ^1 


, = d ♦ d- ; q=1. 

n+1 q 0 ^ 


7. Advance the first and second sums from t to t , 

n n-f I 


v'V , = v*V ♦ f , 

n+1 n n+1 


V^f = V'^f ♦ f , 
n^l n n-*- . 



APPENDIX E 


STARTING ALGORITHM FOR THE GENERAL FORMULATION 


The purpose of a starting procedure for a multistep integra- 
tor is to use the initial conditions (t , y , and the function 

m tn m 

y# y^^^) to approximate the values of f at i nodes where i 
is the order of the integration formula. The Class Il/fixed-mesh gen- 
eral formulation starting algorithm presented here is an iterative pro- 
cedure where the i nodes are at tj^ , k=n , n-1 , . . . ,n-i+1 , t^ is 

t - t 

between t and t . , and t is such that r > 0 , and 

n n-i+1 n h 

n, n-i+1 > 0 . The algorithm can be summarized as follows: 


1. Evaluate f at t : 

m 


f = f(t . y„. y^^^ 
m m m m 


2. Use a Taylor series expansion to obtain the first approximation of 
y , y^^^ and f at the nodes, e.g., 


t. = t + (k-m)h 
K n 


y,, = y^ + (k-m)hy^^^ + (k-m)^h^f 
K n m m 


y^^^ = y_[^^ + (k-m)hf k=n, n-i , . . . ,n-i+1 
K ni ni 


'k = h- »k"> 


where h is the integration steps! ze. 
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f 

i 


13 


3- Form the back differences of f at t to obtain 

n 


r 


.( 1 ). 


V f ; j=1,...,i and form the vector z = (y , y ) 

n p n n 


Approximate the solution and evaluate f at the nodes by using 
the following interpolation formulas: 


fu ^u O) ^2 i ' ,J- 1 _ 

y, = y + (k-m)hy + h ^ a. . V f 

k m m j .j J,k-m n 


(E.1) 


= i 6' _ 


m 


i 

j,k-m n 


j = 1 


(E.2) 


= ^^^k’ ^k’ ^k 


k=n n-i+1 


where 


I 2 2 

a. , = (k-n) a. , - (m-n) a. + (rn-n)(m-k) B. 

j,k-m J,k-n j,m-n Jtm-n 


6. , = (k-n) B- , + (n-m) B^ 

‘^j,k-m J.k-n j,m-n 


and the coefficients a. . and B. are discussed in Appendix 

J • ^ J t K 

A. 


5. Compute the norm of the relative difference between consecutive 

values of (y « ^ 

n n 


2 2 ( S'J* ■ 


where z is the vector or the most recent values of (y , y''^) . 

n n 

Zp is the vector of the previous values of Cy^. 



V 


A 



Zp(j) otherwise 


and q is the number of elements in z 

P 

Finally, if u is not less than some desired tolerance, then the 
iterative procedure is continued by repeating Steps 3i ^ and 5; 
if u is less than the desired tolerance, then the iterative pro- 
cedure is complete and the solution at t^ is computed by setting 
k=n in (E,1) and 

It should be noted that this algorithm may be used for 

any set of nodes relative to t^ . However, the interpolation 

formulas which are used to derive (E.1) and (E.2) are generally 

only valid for (t - t )(t - t . > 0 . In the multistep 

^ n m m n-i^-1 - ^ 

integration packages KSGFS and ABFS, t is chosen so that t 

n m 

is approximately midway between t^ and t^ . Also, the con- 
vergence criterion discussed in Step 5 may be replaced by any 

other appropriate criterion. 



APPENDIX F 


STARTING ALGORITHM FOR THE SECOND-SUM FORMULATION 


The purpose of a starting procedure for a multistep integra- 
tor is to use the initial conditions (t , y , and the function 

in rn !Ti 

r(t, y, to approximate the values of f at the i nodes where 

i is tht order of the integration formula. The Class Il/fixed- 
mesh/second-sum formulation starting algorithm presented here is an 
iterative procedure where the i nodes are at 

t, , k=1, n-1 , . . . ,n-i+1 , t is between t and t . , > 0 . The 
algorithm is summarized as follows: 

1. Evaluate f at t : f = f(t , y , y^^^) 

ra m m m m 

2. Use a Taylor series expansion to obtain the first approximation of 
y , y^^^ , and f at the nodes, e.g.. 


t, = t ♦ (k-m)h 
K n 


2 2 

y = y ♦ (k-m)hy' ' + (k-m) h f 
■'k ■'m 'm n 


( 1 ) ( 1 ) 


y (k-m)hf 


m 


IT. 


''k = '■“k- »k- 


k=n, n-1 , . . . ,n-i+1 


where h is thv* integration stepsize. 
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3. Form the back differences of f at t^ to obtain 




( 1 ). 


V *f^ , j = 1 i and form the vector = (y^, ) 

»l. Approximate the solution and evaluate f at the nodes by using 
the following interpolation formulas: 


( 1 ) .2 


i , 


y , = y + (k-m)hy_ + h 2a.-. 

•'k m m j_.| j-'tX-i 


m n 


(1) (1) u, 5 w ’ f 


fk = y^* yk^^^ • 


where 


(F.1) 


(F.2) 


H 1 . _ = 1 1 , « ” ,m-n 


*j-1,k-m ■ j-1,k-n j-1 ,m-n 


*^j-1,k-m ■ ~ ^j-1,m-n 

and the coefficients a^ ^ and are discussed in Appendix 

B. 


5. Compute the norm of the relative difference between consecutive 


values of (y^^. y„^^) 


2 

D(j) ) 


\( 
J=1 \ 


( 1 ), 


whore z is the vector of the most recent values of (y^* yn ^ • 

, . / (IK 

Zp is the vector of the previous values of (y^, y^ ) 





z (j) otherwise 
P 


and q is the number of elements in 


z 

P 


6, Finally, if u is not less than some desired tolerance, then the 
iterative procedure is continued by repeating Steps 3, ^ and 5; 
if u is less than the desired tolerance, tne iterative procedure 
is complete and the first and second sums are computed at t^ by 


.( 1 ) 


V 'f = . 

n h 


>1 




n n j-i ,m-n n 


It should oe noted that this algorithm may be used for any 

set of nodes relative to t . However, the interpolation formulas 

m 

which are used to derive (F. 1) and (F.2) are generally only valid for 

(t - t )(t - t , ,) > 0 . In the multistep integration packages 
n mm n— i+'i 

SSFSBD and SSFSFE, t is chosen so that t is about midway between 

n m 

and . Also, the convergence criterion in Step 5 may be 

replaced by any other appropriate criterion. 
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